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ABSTRACT 


This report describes the research performed at the Center for Reactive Flow and 
Dynamical Systems in the Laboratory for Computational Physics and Fluid Dynam- 
ics, at the Naval Research Laboratory, in support of the NASA Microgravity Science 
and Applications Program. The primary focus of this research was on investigating 
fundamental questions concerning the propagation and extinction of premixed flames 
in earth gravity and in microgravity environments. 

Our approach was to use detailed time-dependent, multispecies, numerical mod- 
els as tools to simulate flames in different gravity environments. The models include 
a detailed chemical kinetics mechanism consisting of elementary reactions among the 
eight reactive species involved in hydrogen combustion, coupled to algorithms for con- 
vection, thermal conduction, viscosity, molecular and thermal diffusion, and external 
forces. The external force, gravity, can be put in any direction relative to flame prop- 
agation and can have a range of values. 

A combination of one-dimensional and two-dimensional simulations has been used 
to investigate the effects of curvature and dilution on ignition and propagation of flames, 
to help resolve fundamental questions on the existence of flammability limits when there 
are no external losses or buoyancy forces in the system, to understand the mechanism 
leading to cellular instability, and to study the effects of gravity on the transition to 
cellular structure. These studies have shown that a flame in a microgravity environment 
can be extinguished without external losses and that the mechanism leading to cellular 
structure is not preferential diffusion but a thermo-diffusive instability. The simulations 
have also lead to a better understanding of the interactions between buoyancy forces 
and the processes leading to thermo-diffusive instability. When a flame which exhibits a 
cellular structure in a zero-gravity environment is subjected to earth gravity, it evolves 
into either a bubble rising upwards in a tube or a flame oscillating between convex 
and concave curvatures in the case of downward propagation. The effect of gravity is 
minimal on flames in mixtures which are far from flammability limits. 


INTRODUCTION 


In this report we describe the research performed at the Naval Research Labo- 
ratory in support of the NASA Microgravity Science and Applications Program over 
the period December 1985 — December 1988. This work has been performed at the 
Center for Reactive Flow and Dynamical Systems in the Naval Research Laboratory 
by Drs. Gopal Patnaik, K. Kailasanath and Elaine Oran.The emphasis of our research 
has been on investigating fundamental combustion questions concerning the propaga- 
tion and extinction of gas-phase flames in microgravity and earth-gravity environments. 
Our approach to resolving these fundamental questions has been to use detailed time- 
dependent, multispecies numerical models to perform carefully designed computational 
experiments. The basic questions we have addressed, a general description of the nu- 
merical approach, and a summary of the results are described in this introduction. 
More detailed discussions are presented in the next two sections and in the appendices 
to this report. 

It is well known from flammability studies in normal earth gravity that a flame 
propagating upward in a tube propagates in a wider range of mixture stoichiometries 
and dilutions than a flame propagating downward. This means that the apparent 
flammability limits of upward-propagating flames are broader than those of downward- 
propagating flames. Various explanations for this phenomena are based on factors 
such as buoyancy forces, preferential diffusion, flame chemistry, conductive and radia- 
tive heat losses, the aerodynamics of burnt gases, and flame stretch [1-6] . One certainty 
is that gravitational acceleration (buoyancy) is very important and may also be inter- 
acting with and influencing other physical processes. Furthermore, the microgravity 
experiments in the NASA drop tower [7] and Lear-jets [8] have shown that instabilities 
often occur in the propagation of flames near the flammability limits. 

A systematic study which isolates the various processes that could lead to flame 
extinction is needed to gain a better understanding of flammability limits in general 
and more specifically, the role of gravity on flame propagation and extinction. Nu- 
merical simulations, in which the various physical and chemical processes can be inde- 
pendently controlled, can significantly advance our understanding of flame instabilities 
in a microgravity environment and flammability limits on earth and in a microgravity 
environment. In the past three years, we have addressed a number of basic questions 
aimed at resolving some of the important aspects of flame structure and propagation. 
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There axe still many questions left unanswered and, as such, this final report is also a 
progress report on a continuing investigation. 

The basic questions we have adddressed are: the effects of curvature and dilution 
on ignition and propagation of flames, the existence of flammability limits in the absence 
of external losses or buoyancy forces in the system, the mechanisms leading to cellular 
instability, and the effects of gravity on flame instabilities. These studies have shown 
that a flame in a microgravity environment can be extinguished without external losses 
and that the mechanism leading to cellular structure is not preferential diffusion but 
is a thermo-diffusive instability. Furthermore, the simulations have also lead to a 
better understanding of the interactions between gravity induced flame dynamics and 
the intrinsic processes involved in flame propagation. When a flame which exhibits a 
cellular structure in a zero-gravity environment is subjected to earth gravity, it evolves 
into either a bubble rising upwards in a tube or a flame oscillating between convex and 
concave curvatures in the case of downward propagation. The simulations have also 
shown that the effect of gravity is minimal on flames in mixtures which are far from 
flammability limits. 

In the research described here, we have used both one-dimensional and two- 
dimensional numerical simulations to systematically isolate and evaluate the impor- 
tance of various processes that might be controlling the dynamics of flames, particularly 
near the flammability limits, and to evaluate the relative importance of these processes 
in normal gravity and microgravity. Both the numerical models axe time-dependent 
and solve the multispecies coupled partial differential reactive-flow equations. These 
models include a detailed chemical kinetics mechanism coupled to algorithms for con- 
vection, thermal conduction, viscosity, molecular diffusion, thermal diffusion, and ex- 
ternal forces. The external force, gravity, can be in any direction relative to flame 
propagation and can have a range of values. Energy sources and sinks may also be 
prescribed as a function of time. The one-dimensional model was developed before this 
NASA project began, but the two-dimensional model was developed during the past 
two years. The development of the two-dimensional model was supported primarily 
by ONR/NRL funding, and it has been applied to study the evolution of multidimen- 
sional flame structure in support of the NASA Microgravity Science and Applications 
Program. The models are described in greater detail in the next section. 
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COMPUTATIONAL TOOLS 


Two computational models developed at NRL, FLIC2D and FLAME1D, have been 
used to study the propagation and extinction of premixed gaseous flames. FLIC2D has 
been used to study the tendency of hydrogen-air mixtures to show cellular instability 
[9] and more recently to model diffusion flames [10]. Many submodels in FLIC2D axe 
based on FLAME1D, a one-dimensional, time-dependent Lagrangian model which has 
been tested extensively and applied to various studies of flame phenomena including 
calculations of burning velocities [11,12], minimum ignition energies and quenching 
distances [12,13], effects of curvature and dilution [14], and flammability limits [15]. 
For completeness, both models are described below. 

FLAME1D 

This is a time-dependent model [12] which solves the compressible, conservation 
equations for mass, momentum and energy [16,17] in one spatial dimension. Since it 
was developed specifically to study the initiation, propagation, and quenching of lami- 
nar flames, it consists of algorithms for modelling convection, detailed chemical kinetics 
and energy release, thermal conduction, molecular and thermal diffusion of the individ- 
ual species, and various energy deposition mechanisms. The model has a modular form 
and the algorithms representing the various chemical and physical processes are com- 
bined using an asymptotic timestep-split approach in which the individual processes 
are integrated separately and then coupled together [17,18]. The model also permits a 
variety of initial and boundary conditions. 

The convective transport terms in the equations are solved by the algorithm AD- 
INC, a Lagrangian convection algorithm which solves implicitly for the pressures [19]. 
Since ADINC communicates compression and expansion across the system implicitly,, 
it overcomes the Courant time-step limit. Since it is Lagrangian, it can maintain steep 
gradients computationally without numerical diffusion for a long period of time. This 
is important in flame calculations where the diffusive transport of material and energy 
can govern the system evolution and therefore must be calculated accurately. An adap- 
tive regridding algorithm has been developed to add and delete computational cells. 
For the flame calculation, cells are added in the region ahead of the flame so that the 
front always propagates into a finely zoned region. Cells are removed behind the flame 
where the gradients are very shallow. 
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The chemical interactions are described by a set of nonlinear coupled ordinary dif- 
ferential equations. This set of equations may be stiff when there are large differences 
in the time constants associated with different chemical reactions, and are invariably 
stiff for combustion problems. These equations are solved using VSAIM, a fully vec- 
torized version of the selected asymptotic integration method CHEMEQ [20,21]. In 
VSAIM, the stiff equations are identified and solved using a very stable asymptotic 
method while the remaining equations are solved using a standard classical method. 

The diffusive transport processes considered in this model are molecular diffusion, 
thermal conduction and thermal diffusion. These processes are crucial to the description 
of flame phenomena since they are the mechanisms by which heat and reactive species 
a re transported ahead to the unburned gas. An iterative algorithm, DFLUX, is used to 
obtain the diffusion velocities without the cost of performing matrix inversions [17,22]. 
This method has been vectorized and is substantially faster than matrix inversions 
when four or more species are involved. 

Problems can be set up in either planar, cylindrical, or spherical coordinates as 
well as a variable power series coordinates which can model one-dimensional nozzle-like 
geometries. The model can also be configured with either an open or closed boundary 
at one end. The open boundary simulates an unconfined system. 

FLIC2D 

FLIC2D is a time- dependent, Eulerian, implicit, compressible, two-dimensional 
flame model. In FLAME1D, a Lagrangian method was chosen for convective transport 
because eliminating the advection term in effect means eliminating numerical diffusion 
from the calculation. Extending, this Lagrangian approach to multidimensions would 
be extremely difficult and expensive, requiring many years of algorithm development 
and testing. Therefore, we felt that we would have more success in multidimensions 
with an Eulerian method. However, most Eulerian methods are either more numer- 
ically diffusive than what is acceptable in a flame model, or they are explicit and 
hence extremely inefficient at the very low velocities associated with laminar flames. 
To circumvent these numerical problems, we developed BIC-FCT, the Barely Implicit 
Correction to Flux- Corrected Transport [23,24]. BIC-FCT combines an explicit high- 
order, nonlinear FCT method [25,26] with an implicit correction process. This combi- 
nation maintains high-order accuracy and yet removes the timestep limit imposed by 
the speed of sound. By using FCT for the explicit step, BIC-FCT is accurate enough 
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to compute with sharp gradients without overshoots and undershoots. Thus spurious 
numerical oscillations that would lead to unphysical chemical reactions do not occur. 
The development of this new algorithm has made it possible to model multidimensional 
flames in detail. However, because there is always some residual numerical diffusion, 
even in high-order Eulerian algorithms, we need to monitor calculations to ensure that 
numerical diffusion never becomes larger than the physical diffusion processes we need 
to resolve. 

In generalizing the flame model to two-dimensions, simplifications were made in 
the diffusive transport calculations in order to reduce the cost of computations. In 
FLIC2D, thermal conductivity of the individual species is modeled by a polynomial 
fit in temperature to existing experimental data. Individual conductivities are then 
averaged using a mixture rule [12,27] to get the thermal conductivity coefficient of the 
gas mixture. A similar process is used to obtain the mixture viscosity from individual 
viscosities. Heat and momentum diffusion are then calculated explicitly using these 
coefficients. 

Mass diffusion also plays a major role in determining the properties of laminar 
flames. Binary mass-diffusion coefficients are represented by an exponential fit to exper- 
imental data, and the individual species-diffusion coefficients are obtained by applying 
mixture rules [12] .The individual species-diffusion velocities are determined explicitly 
by applying Fick’s law followed by a correction procedure to ensure zero net flux [27], 
This procedure is equivalent to using the iterative algorithm DFLUX [22] (used in 
FLAME1D) to second order. 

The chemical reaction-rate equations are modelled and solved as before using a 
vectorized version of CHEMEQ, an integrator for stiff ordinary differential equations 
[21]. Because of the complexity of the reaction scheme and the large number of com- 
putational cells in a two-dimensional calculation, the solution of the chemical rate 
equations takes a large fraction of the total computational time. A special version of 
CHEMEQ called TBA was developed to exploit the special hardware features of the 
CRAY X-MP. 

As in the one-dimensional code, all of the chemical and physical processes are 
solved sequentially and then coupled asymptotically by timestep splitting [17,18]. This 
modular approach greatly simplifies the model and makes it easier to test and change 
the model. Individual modules were tested against known analytic and other previously 
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verified numerical solutions. One-dimensional predictions of the complete model were 
compared to those from the Lagrangian model FLAME1D which has been benchmarked 
extensively against theory and experiment. 

In summary, all of the chemical and diffusive-transport algorithms are essentially 
the same in FLIC2D as in FLAME1D. Thus given the input data concerning the species, 
their transport coefficients, and the reactions among the species, either the one- or 
two-dimensional model can be run. The major differences between the two models 
jure in the costs of running them and the algorithm chosen for convective transport. 
In general, we use the results of one-dimensional calculations as initial conditions for 
the two-dimensional calculations in order to reduce the cost of computations. The 
two-dimensional flame model is discussed in greater detail in Appendix A. 


SUMMARY OF RESEARCH 


The thrust of the work during the past three years was two-fold. The first part 
was to study the propagation and extinction of flames using FLAME1D, the time- 
dependent, one-dimensional flame model. The second was to develop and use a mul- 
tidimensional flame model, FLIC2D, to study mechanisms leading to cellular flame 
instability and the effects of gravity on the multidimensional structure of flames. Both 
models use a detailed set of chemical kinetics of all the species involved in hydrogen 
combustion as well as multi-species difffusion and thermal conduction. 

Effects of Curvature and Dilution on Flame Propagation 

We studied the effects of curvature and dilution on zero-gravity flame propagation 
in hydrogen-oxygen- nitrogen mixtures using the time-dependent, one-dimensional, La- 
grangian flame model, FLAME1D. We performed a systematic study in which we varied 
both dilution and geometry. First, we considered the propagation of a spherical flame 
in a stoichiometric hydrogen-air mixture (2:1:4). We then increased the amount of 
nitrogen dilution and simulated spherical flames in hydrogen-oxygen-nitrogen mixtures 
in the ratios 2:1:7, 2:1:10, 2:1:11, 2:1:13 and 2:1:15 . We finally studied the propagation 
of cylindrical and planar flames in the same mixtures. 

We found that in a stoichiometric hydrogen-air mixture, the velocity of a spherical 
flame first decreases with increasing radii until it reaches a minimum value, and then it 
increases. For large radii, the burning velocity approaches the planar burning velocity. 
With increasing dilution, the same trends were observed. Furthermore, the burning 
velocity of both spherical and planar flames decreases with increasing dilution. 

These observations lead to the conclusion that a flame can be extinguished 
with less dilution in one geometry than in another. Specifically, we showed that a 
2:1 :13/hydrogen-oxy gen- nitrogen mixture does not support a self-sustained spherical 
flame, although it does support cylindrical and planar flames. These observations 
have been explained on the basis of preferential diffusion, flame stretch, and chemical 
kinetics. These observations are also in qualitative agreement with experimental ob- 
servations in both zero and one gravity. Detailed comparisons with experimental data 
are not possible because the actual flames near this extinction limit are unstable and 
multidimensional. 
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We are currently extending this work to fuel- lean and fuel- rich mixtures. Pre- 
liminary results indicate that the same trends are observed in fuel-rich mixtures as in 
stoichiometric mixtures discussed above. However, in fuel- lean mixtures, we expect the 
flame velocity to first increase with increasing radii and then to decrease to the planar 
burning velocity for large radii. This might result in mixtures which burn for some time 
and then extinguish for large radii where the burning velocity drops. Such behavior is 
characteristic of the “self-extinguishing” flames observed in the NASA drop tower [7]. 
In the actual experiment, heat loss might play an important role [28]. Furthermore, 
lean hydrogen-air mixtures are prone to cellular instabilities. Therefore, detailed study 
of flame propagation in such mixtures must be done with a multidimensional model. 

A paper describing our work in this area in greater detail is enclosed as Appendix B. 

Flame Propagation and Extinction in Fuel-Rich Mixtures 

Microgravity flames in fuel-rich hydrogen-air mixtures appear to be stable and 
nearly one-dimensional. Therefore, we can study the extinction process and flamma- 
bility limits by simulating the propagation of planar flames in fuel-rich hydrogen-air 
mixtures at zero gravity. We found that the burning velocity and flame temperature 
monotonically decreased with increasing fuel concentration. With further increase in 
the fuel-concentration, the temperature of the burning mixture does not seem to be 
high enough to sustain a steadily propagating flame. For temperatures of about 950 K, 
the endothermic radical-producing reactions take up more energy than what is pro- 
duced by the exothermic reactions and hence a self-sustained flame is not observed. In 
a practical system, even minimal heat losses would ensure that an incipient flame is 
extinguished in such a fuel-rich mixture. The exact limit composition and temperature 
would vary from one experimented set-up to another. This work is described in greater 
detail in Appendix C. 

Simulations of Cellular Flames 

After testing the FLIC2D code, we used it to simulate unsteady nonplanar flames 
and the mechanisms leading to cellular structure. The first calculation modelled a flame 
in a fuel-lean mixture of hydrogen-oxygen-nitrogen in the ratio of 1.5:1:10. A multidi- 
mensional flame is observed in this mixture in the experiments of Mitani and Williams 
[29]. The numerical simulations showed that in this mixture, a one-dimensional planar 
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flame is unstable because perturbing the flame leads to a multidimensional structure 
resembling that of a cellular flame. 

A similar calculation was performed for a fuel-rich mixture of hydrogen-oxygen- 
nitrogen in the ratio, 3:1:16. Flames in this mixture are stable in the experiments 
of Mitani and Williams [29], and also in our calculations. In this case, an initial 
perturbation is quickly damped, restoring the flame profile to its one-dimensional shape. 

We then performed additional simulations to determine the mechanism controlling 
the instability. There are two possible mechanisms that can lead to cellular instabil- 
ity. One mechanism is the preferential diffusion of one reactant over another due to 
its higher diffusivity. For example, in a hydrogen-oxygen-nitrogen mixture, hydrogen 
diffuses much faster than oxygen or nitrogen. Another mechanism is the differences 
in the rate at which mass and heat diffuse (diffusional-thermal theory). In the third 
calculation, the mass diffusivity of hydrogen was set equal to oxygen. This eliminates 
the preferential diffusion effect, and means that the light fuel will not move fast com- 
pared to the other reactant. If the mechanism of instability were preferential diffusion, 
the lean flame would be stable. This is in fact the result: with this change of relative 
diffusivities, the lean flame was stable in the calculation. 

This result supports both the preferential diffusion theory and the diffusional- 
thermal theory. A further calculation was performed in which the mass diffusivity 
of oxygen was set equal to that of hydrogen. In this case, the flame was unstable, 
which disagrees with the prediction of preferential diffusion. However, the diffusional- 
thermal theory does predict that this flame is unstable. Thus for these hydrogen flames 
studied, the results agree with the predictions of the diffusional-thermal theory and not 
the preferential diffusion theory. These simulations provide an example of the use of 
numerical simulations in identifying the controlling instability when more than one 
mechanism is plausible. A paper describing detailed numerical simulations of cellular 
flames is attached to this report as appendix D. 

Effects of Gravity on Cellular Flames 

We have studied the effects of gravity on flame structure by comparing simulations 
of zero-gravity flames to upward- and downward-propagating flames. These simulations 
show that the effects of gravity become greater as the lean flammability limit is ap- 
proached. For example, in a 1.5:l:10/hydrogen-oxygen-nitrogen mixture, gravity plays 
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only a secondary role in determining the multidimensional structure of the flame. How- 
ever, in a 1:1:10 mixture, an upward- propagating flame is highly curved and evolves 
into a bubble that rises upwards in the tube. The zero-gravity flame shows a cel- 
lular structure. The structure of the downward-propagating flame oscillates in time 
exhibiting both concave and convex curvatures towards the unburnt mixture. These 
observations have been explained on the basis of an interaction between the processes 
leading to Rayleigh- Taylor instability and the cellular instability. The simulations also 
suggest that cellular instability grows more rapidly than Rayleigh- Taylor instability. A 
paper describing this work in greater detail is included as appendix E. 

A list of papers and presentations related to the work described above are given 
below. 

Journal Articles and Meeting Proceedings 

Effects of Curvature and Dilution on Unsteady Flame Propagation, K. Kailasanath and 
E.S. Oran, in Dynamics of Reactive Systems, I. Flames and Configurations, Prog. 
Aero. Astro. 105, AIAA, 1986. 

A Barely Implicit Correction for Flux-Corrected Transport, G. Patnaik, R.H. Guirguis, 

J. P. Boris, and E.S. Oran, J. Comput. Phys., 71, 1-20 (1987). 

An Implicit Flux- Corrected Transport Scheme for Low Speed Flow, G. Patnaik, J.P. 
Boris, R.H. Guirguis, and E.S. Oran, AIAA 25th Aerospace Sciences Meeting, 
Paper No. AIAA-87-0482, AIAA, NY. 

Time- Dependent Simulations of Laminar Flames in Hydrogen-Air Mixtures, K. 
Kailasanath and E.S. Oran, in Complex Chemical Reaction Systems, Eds., J. 
Warnatz and W. Jager , Springer- Verlag, Heidelberg, 1987. 

Numerical Modelling of Unsteady Laminar Flames, K. Kailasanath, E.S. Oran and J.P. 
Boris, to appear as Chapter 3 in Physics, Chemistry and the Modelling of Flames 
by R. Fristrom, Johns Hopkins University Press, 1988. 

Simulations of Flames near the Rich Flammability Limit of Hydrogen- Air Mixtures, 

K. Kailasanath and E.S. Oran, in preparation for submission to Combustion and 
Flame, 1989. 

Detailed Numerical Simulations of Cellular Flames , G. Patnaik, K. Kailasanath, K. J. 
Laskey, and E. S. Oran, to be published in the Proceedings of the 22th Symposium 
(International) on Combustion, Seattle, WA., August, 1988. 


10 



FLIC: A Detailed Two-dimensional, Flame Model, G. Patnaik, K. J. Laskey, K. 
Kailasanath, E.S. Oran and T.A. Brun, to be published as Naval Research Labo- 
ratory Memorandum Report, 1989. 

Meeting Presentations 

Effects of Geometry and Diffusive Transport on Premixed, Laminar Flame Propagation, 
K. Kailasanath, E.S. Oran and G. Patnaik, 1985 Fall Technical Meeting of the 
Eastern Section of the Combustion Institute, Philadelphia, November, 1985. 

Effects of Fluctuations on Chemically Reacting Systems and the Generation of Tur- 
bulence, E. S. Oran, invited talk presented at the Workshop on Flame Structure, 
Novosibirsk, USSR, July, 1986. 

Ignition Studies of a Dilute, Stoichiometric Mixture of Hydrogen, Oxygen and Nitrogen, 
S. Kurban Ghafil, C. Brochet, K. Kailasanath, and E. Oran, Poster Presentation at 
the 21th Symposium (International) on Combustion , Munich, Germany, August, 
1986. 

Time-Dependent Simulations of Laminar Flames in Hydrogen-Air Mixtures, K. 
Kailasanath and E.S. Oran, Second Workshop on Modelling of Complex Reac- 
tion Systems, Heidelberg, Germany, August, 1986. 

Application of an Implicit Flux- Corrected Transport Scheme to Combustion Modeling, 
G. Patnaik, K. Kailasanath, and E. Oran, 1986 Fall Technical Meeting of the 
Eastern Section of the Combustion Institute, San Juan, PR, December, 1986. 

An Implicit Flux- Corrected Transport Scheme for Low Speed Flow, G. Patnaik, J.P. 
Boris, R.H. Guirguis, and E.S. Oran, AIAA Paper 87-0482, 25th Aerospace Sci- 
ences Meeting, AIAA, Reno, January, 1987. 

Some Challenges in the Numerical Simulation of Laminar and Turbulent Flames, Elaine 
S. Oran, Princeton University, March, 1987 

Numerical Calculations of Combustion Flows, Elaine S. Oran, University of California, 
San Diego, November, 1987. 

Detailed Numerical Simulations of Lean Flame Structure, G. Patnaik, K. Kailasanath, 
K. J. Laskey, and E. S. Oran, 1987 Fall Technical Meeting of the Eastern Section 
of the Combustion Institute, Gaithersburg, MD, November, 1987. 

Numerical Simulations of Cellular Structure in Hydrogen Flames, G. Patnaik, K. 
Kailasanath, K. J. Laskey, and E. S. Oran, Presented at the 3rd Workshop on 
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Numerical Methods in Laminar Flame Propagation, Leeds, England, April, 1988. 

Numerical Simulation of Complex Fluid Systems, Elaine S. Oran, George Washington 
University, Washington, DC, April, 1988. 

Future Trends in Numerical Simulations of Reacting Flows, Elaine S. Oran, 12th In- 
ternational Association for Mathematics and Computers in Simulation (IMACS) 
World Congress, Paxis, France, July, 1988. 

Detailed Numerical Simulations of Cellular Flames , G. Patnaik, K. Kailasanath, K. 
J. Laskey, and E. S. Oran, presented at the 22th Symposium (International) on 
Combustion, Seattle, WA., August, 1988. 

Effect of Gravity on Multi-Dimensional Laminar Premixed Flame Structure, K. 
Kailasanath, G. Patnaik and E. S. Oran, presented at the 39th International As- 
tronautical Congress , Bangalore, India, October, 1988. 

New Directions in Computing Reacting Flows, Elaine S. Oran, Symposium on Ad- 
vances and Trends in Computational Structural Mechanics and Fluid Dynamics, 
Washington, DC, October, 1988. 

Numerical Simulations of Unsteady, Unstable Flames, E.S. Oran, G. Patnaik, K. 
Kailasanath, and K. Laskey, American Physical Society, Division of Fluid Dy- 
namics, November, 1988. 

The Current Status of CFD in Combustion Modeling, E. S. Oran, invited talk to 
be presented at the 1988 Fall Technical Meeting of the Eastern Section of the 
Combustion Institute, Clearwater, FL, December, 1988. 
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FUTURE DIRECTIONS 


The calculations discussed above have lead to a better understanding of the struc- 
ture and propagation of flames. We now understand the reasons for some of the ob- 
served differences in the structure of upward and downward propagating flames. How- 
ever, we need to do further simulations to completely describe upward and downward 
propagating flames in the standard flammability limit tube and to determine flamma- 
bility limits. An important factor that has not been considered in the above simulations 
is the effect of losses due to heat conduction to walls and due to radiation. Furthermore 
flames near the flammability limits are actually three-dimensional. In principle, FLIC 
can be extended to perform calculations in a three-dimensional geometry. Currently, a 
major limiting factor is the cost of such computations. 

The solution of the chemical rate equations takes a significant portion (60 - 70 %) 
of the computer time because of the large number of species and the large number 
of rate equations involved in a mechanism consisting of elementary reactions. The 
problem becomes even more expensive in hydrocarbon fuels, starting with the simplest, 
methane. One approach to reducing the cost of multidimensional flame simulations is to 
use a greatly simplified reaction mechanism. However, a satisfactory simplified reaction 
mechanism does not exist currently, even for methane. Although, we have reduced a 108 
reaction-rate mechanism for methane to about 50 rates [30], this is still too expensive to 
use in a multidimensional flame calculation. Currently there are many ongoing efforts 
to find ways of easing the burden of integrating full chemical reaction-rate mechanisms 
[31,32]. Using such simplified approaches in FLIC2D is straightforward and will enable 
us to simulate multidimensional flames in hydrocarbon mixtures. 

Another approach to reducing the cost of chemistry calculations is to develop 
faster chemical integration algorithms. We are continually finding ways to make the 
algorithms we use faster. In addition, there are some new approaches to sorting stiff 
equations [33] which are promising. However, it is unlikely that reducing the computer 
time required by using faster algorithms alone will be sufficient to allow us to perform 
calculations of complex hydrocarbon flames in the near future. 

The flame calculations performed to date with FLIC2D are in a two-dimensional 
cartesian coordinate system. With this version of the code we have been able to simulate 
the cellular structure of flames in a channel. Extensions to an axisymmetric geometry 
are straightforward. The axisymmetric version of the code will be used to simulate flame 
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propagation in cylindrical tubes. This will enable us to simulate new modes of flame 
instability. Calculations in an axisymmetric geometry using simplified and appropriate 
chemistry need to be done before full three-dimensional simulations of flames in tubes 
are performed. 
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Abstract 


This report describes FLIC, a detailed, two-dimensional model for flames. In or- 
der to describe a flame in enough detail to simulate its initiation, propagation, and 
extinction, FLIC combines algorithms for subsonic convective transport with buoy- 
ancy, detailed chemical reaction processes, and diffusive transport processes such 
as molecular diffusion, thermal conduction, and viscosity. Several new numerical 
techniques had to be developed specifically for multi-dimensional flame modelling, 
and these are highlighted in this report. The numerical implementation of the flame 
model is described here, and a typical application is provided for illustration. FLIC 
is currently the state-of-the-art in detailed two-dimensional, time dependent flame 
models. 
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1. Introduction 


This report describes the new computer code FLame with Implicit Convection 
(FLIC), a two-dimensional time-dependent program developed specifically to com- 
pute and study the behavior of flames and other subsonic chemically reactive flows. 
In order to describe a flame in enough detail to simulate its initiation, propagation, 
and extinction, FLIC combines algorithms for subsonic convective transport with 
buoyancy, detailed chemical reaction processes, and diffusive transport processes 
such as molecular diffusion, thermal conduction, and viscosity. Currently, we have 
not included algorithms for radiation transport or thermal diffusion, although these 
important processes can, in principle, be added in the same modular fashion as 
those physical processes that have been included. 


FLIC solves the reactive-flow conservation equations for density, p, momentum, 
pV, energy, E, and number densities of the individual species, n k , k = 1, ...,n 4p , 
according to: 


| + V.(,V) 

^r +v -('™) 

f + v.(eV) 

IT + 


0, 

(1.1) 

-VP + F-VxpVxV + V (|/iV • Vj , 

(1.2) 

-V • (PV) + V • (acVT) - 


n#p n r 

'£V-(n l h l V t )+£<?., 

*=1 r=l 

(1.3) 

-V • ( n k V k ) + w k . 

(1.4) 


Here V is the fluid velocity, P is the pressure, p is the coefficient of shear viscosity, 

F is a body force, « is the thermal conductivity of the mixture of gases, h k is 

— ♦ 

the enthalpy of species k, V k is the diffusion velocity of species k, Q r is the heat 
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released from reaction r, and w k is production of species k by chemical reaction. 
These equations are solved assuming that the individual species are ideal gases 
obeying the thermal equation of state, 

Pk=n k kT , (1.5) 

and that the differential relation between internal energy u and pressure P is given 
by 

Su = — ( 1 . 6 ) 

7-1 

where 7, the ratio of specific heats of the mixture, is a function of its temperature 
and composition. 

FLIC solves these equations in two dimensions for low-velocity compressible 
flow such that the Mach number is less than 0.1. Depending on specified initial and 
boundary conditions, the flames modeled can be either premixed or diffusion flames 
in either planar or axisymmetric geometries. 

To date, FLIC has been applied to the study of the cellular instability near 
the flammability limits of a premixed flame in a gas containing hydrogen, oxygen, 
and nitrogen flame [1]. Cellular flames are formed when the mass diffusion of the 
deficient reactant overwhelms the stabilizing influence of heat conduction [2,3]. For 
hydrogen flames, this behavior is seen close to the extinction limit [4] when the flame 
is thick and temperatures are low. The fluid velocity is usually low, typically less 
than 20 cm/s in the burned region. Radiation effects are not important because 
the flame is not luminous and absorbing species, such as CO2 found in hydrocarbon 
flames, are absent. A brief summary of this work is given in section 6. 


1.1 Algorithm Development and Implementation 


Producing a code that describes low-speed flames required the development of sev- 
eral new numerical methods as well as finding new ways to implement existing 
algorithms. For example, a new multidimensional, low-speed convection algorithm 
had to be developed. An important requirement of any convection algorithm is 
that the numerical diffusion not be larger than the physical diffusion processes that 
must be resolved. The FCT method [5,6] meets this criterion, but it is inherently 
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an explicit method; so that the small timesteps it requires makes it inefficient for 
low-speed flows. The BIC-FCT method [7] was developed specifically to combine 
the accuracy of FCT with the efficiency of an implicit method for low-speed flows. 
BIC-FCT allows timesteps up to a hundred times larger than FCT and yet the cal- 
culation time of one BIC-FCT timestep is approximately equal to the calculation 
time of one timestep in the standard FCT module. Section 2 describes BIC-FCT 
in detail. 

To solve the detailed combustion equations for hydrogen-oxygen chemistry in- 
volves solving coupled nonlinear ordinary differential equations for eight species and 
48 chemical reactions representing the conversion of chemical species and chemical 
energy release into the system. This is the most expensive part of a reactive flow 
calculation because it requires integrating the set of equations at each computa- 
tional cell for each timestep. The characteristic times of these differential equations 
vary by orders of magnitude, resulting in a set of very “stiff” equations. Then, 
because the cost of the calculation is approximately linear with the number of 
computational cells, the computational cost can be extreme in multidimensional 
computations. FLIC handles the cost of integrating ODEs in two ways: one, by 
not integrating the chemical reaction equations where there is nothing or essentially 
nothing happening, the other is by optimizing the integration procedure. We are 
using the CHEMEQ [8,9] method to solve stiff sets of ordinary differential equations, 
but in the TBA implementation that is fully optimized for the CRAY X-MP com- 
puter. TBA, described in section 3, allows speeds of up to 50 percent over VSAIM, 
the multidimensional implementation of CHEMEQ used to date. 

Thermo-physical properties of the individual species and the mixture are re- 
quired throughout the computation. These properties are modelled with high-order 
curve fits to values derived from more accurate calculations. The individual proper- 
ties are combined where needed to obtain mixture properties through mixing rules. 
This simplified method is highy efficient yet sufficiently accurate. Section 4 describes 
the numerical solution of the diffusive transport processes. 

The submodels representing the various physical processes are in independent 
modules that are coupled together. Several modifications have been made to the 
usual timestep splitting method in order to increase the stability limits and improve 
the efficiency of FLIC. Section 5 describes the details of these improvements. 
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1.2 Comparison with FLAME1D 


FLAME1D is a one-dimensional program used extensively to study the properties 
of the ignition and extinction of hydrogen flames [10]. It is useful to compare FLIC 
to FLAME1D because some FLAME1D algorithms are not obviously extendable to 
multidimensions and others are simply too expensive. The more important differ- 
ences between these codes are described here. 

1. FLIC uses an Eulerian representation of the convective transport instead of 
a Lagrangian representation. Specifically, ADINC, the one-dimensional La- 
grangian algorithm [11] used in FLAME1D, is replaced by BIC-FCT [7], a very 
accurate multidimensional implicit Eulerian algorithm. A Lagrangian formu- 
lation, though preferable, is exceedingly difficult in multidimensions. However, 
unlike ADINC, BIC-FCT can be readily used to describe two-dimensional or 
three-dimensional flows. 

2. Although the basic CHEMEQ algorithm is used in both codes, the VSAIM 
implementation used in FLAME1D was replaced by the TBA implementation. 
This algorithm is optimized for the CRAY X-MP, and can be retrofitted into 
CRAY versions of FLAME1D. 

3. Another expensive part of the flame calculation is determining the amount 
that individual species diffuse. In FLAME1D, we used an iterative matrix 
expansion algorithm [12] that produces the diffusion flux of a species to arbi- 
trary order, although we generally used it only up to second order. In FLIC, 
we use a simpler technique which first evaluates a Fickian flux and then makes 
a correction. This approach is equivalent to the first order of the matrix ex- 
pansion, cannot be made higher order, but it is computationally less intensive 
and certainly adequate for the flame problems treated to date. 

4. In FLAME1D, thermal conduction is computed directly from expressions for 
the individual thermal conductivities of the species derived from molecular 
theory. In FLIC , we use curve fits and mixture rules which has been bench- 
marked against the more exact calculations used in FLAME1D. 

5. Whereas FLAME1D did not include algorithms for either viscosity or gravity, 
both of these effects are included in FLIC. 
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6. At the moment, FLAME1D has a more flexible gridding algorithm. Because 
the basic convection algorithm in FLAME1D was Lagrangian, relatively non- 
diffusive cell splitting and merging routines are used to refine or coarsen the 
grid. The result is a rather general gridding capability. The general approach 
to adaptive gridding must be different in a multidimensional Eulerian code, 
and this is now being developed for FLIC. 


5 


A-9 


2. Convection 


In this chapter, we describe the fluid convection algorithm, BIC-FCT, how it is 
used in FLIC , and how the gravitational acceleration term is included. In detailed 
flame simulations that must resolve the individual species diffusion, the numerical 
diffusion that results from solving the convection equations numerically must be 
small enough so that we can resolve the physical diffusion processes. For high- 
speed flows, the Eulerian explicit monotone methods such as FCT [13] or most 
of the TVD [14] methods achieve this goal, and some even allow variable-order 
accuracy. Unfortunately, the timestep required by explicit methods must be small 
enough to resolve the sound waves in the system, otherwise the numerical method 
is unstable. There is little or no penalty paid for this small timestep in high-speed 
flow in which the physical phenomena evolve fast, but for low-speed flows, explicit 
methods are very inefficient and expensive. For example, resolving a microsecond 
of physical time with a timestep of 10~ 8 s requires 100 timesteps, but resolving one 
second requires 10 8 timesteps. For many low-speed flows this temporal accuracy 
is unnecessary; we need only to be able to resolve a millisecond of physical time 
with 100 timesteps. For this reason, implicit methods that allow large timesteps 
are usually used for calculations of low-speed flows. 

One often-considered approach to eliminating numerical diffusion is to use La- 
grangian methods, in which diffusion is totally absent by definition. We have found 
that this approach works well in one dimension, but there are a number of seri- 
ous problems in multi-dimensions [15]. In complex flows, the multidimensional La- 
grangian grid becomes distorted to the point where nearest neighbors are no longer 
connected by grid lines, a situation that leads to extremely inaccurate calculations. 
Eventually grid lines can even cross, which makes the solution unstable. Such prob- 
lems are often avoided by a regridding procedure that actually adds diffusion to 
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the solution, or by using dynamically restructuring grid such as in SPLISH [16], 
which adds complexity. Except in one dimension, it has not yet been shown that 
Lagrangian methods are to be preferred because of the geometric complexities that 
arise. 

Most common methods for solving convection problems use algorithms that 
produce ripples near steep gradients such as in a flame or shock front. The first high- 
order, nonlinear, monotone algorithm, Flux-Corrected Transport (FCT)[13], was 
designed to prevent these ripples by maintaining local positivity near steep gradients 
while keeping a high order of accuracy elsewhere. Other nonlinear methods have 
been reviewed by Woodward and Collela[14]. Although these methods are explicit, 
there are recent reports on implicit, nonlinear methods [17,18]. A major problem 
with applying these implicit methods to low-speed flows is that they are expensive 
even though they can be very accurate. 

The Barely Implicit Correction to Flux- Corrected Transport, BIC-FCT[7], was 
designed to overcome the problem of numerical diffusion in low-velocity implicit 
methods. BIC-FCT combines an explicit high-order, nonlinear FCT method [5,6] 
with an implicit correction process. This method removes the timestep limit im- 
posed by the speed of sound on explicit methods, retains the accuracy required to 
resolve the detailed features of the flow, and keeps the computational cost as low 
as possible. 


2.1 The BIC-FCT Algorithm 

BIC-FCT is based on am approach suggested by Casulli and Greenspan [19], who 
showed that it is not necessary to treat all of the terms in the gas-dynamic equations 
implicitly to be able to use longer timesteps than those dictated by explicit stability 
limits. Only those explicit terms which force a timestep limit due to the sound speed 
have to be treated implicitly. In a pure convection problem, the timestep is still 
limited by the fluid velocity, but for the low Mach number flow in flames, this results 
in a hundredfold increase in the timestep. The term “Barely Implicit Correction” 
emphasizes that only the minimal number of terms in the conservation equations 
are treated implicitly. 
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BIC-FCT solves the convective portion of the Navier-Stokes equations: 


| + v -(,v) = 0, 

(2.1) 

+ V • (pVV) VP + F, 

(2.2) 

^ + v(bv) = -V ■ (pv) + 5 , 

(2.3) 

ir+ v -M) = 0. 

(2.4) 


where F is a body force, usually gravity, and S is used to couple in the contributions 
to the change in energy due to other processes as discussed in chapter 5. In addition, 
note that convection of species is done at the same time. 


BIC-FCT takes a predictor-corrector approach. The first explicit step uses FCT 
and a large timestep governed by a CFL condition on the fluid velocity, 


P'~P° 

At 

p'V' - p°V° 
At 

E' - E° 
At 


-V • (p°V°) , 

- V • ( P°V°V °) - VP°, 

-V • (E° + P°) [u>V f + (1 - u)V°\ + 5, 


(2.5) 

( 2 . 6 ) 
(2.7) 


where the implicitness parameter 0.5 < u> < 1.0. This produces the intermediate 
values denoted by primes. 


The second step is an implicit correction requiring the solution of one elliptic 
equation for the pressure correction, 6P = u(P n — P°) : 


6P 

(7 — l)u;At 


— u>A £V 


' E° + P° 

, ft 


V6P\ = 


E' - E° 
At 


p'V' 2 - p°V 2 

2At 


( 2 . 8 ) 


The elliptic equation (Eq. (2.8)) is solved by the multigrid method, MGRID [20]. 
This method is O(N) in both number of computations as well as storage. It is 
vectorized and extremely efficient on the CRAY X-MP. MGRID requires that the 
number of grid points in each direction be factorizible by a laxge power of two. 
While this has not proven to be very restrictive for FLIC, it has been a problem 
in other applications. MGRID also has only a limited set of boundary conditions. 
Fortunately, the Neumann boundary conditions used in solving Eq. (2.8) has been 
implemented. 


8 


A- 12 



The final step is the correction of the provisional values for momentum, energy, 
and pressure: 


p n v n 

= -A tVSP + p'V’, 

(2.9) 

E n 

SP 

— , , x + F? , 

( 7“ IV 

(2.10) 

pn 

= 6P + P 0 . 

(2.11) 


Because BIC-FCT uses FCT for the explicit step, it has the high-order monotone 
properties that accurately treat sharp gradients. This accuracy combined with the 
savings that result from removing the sound-speed limitation on the timestep makes 
BIC-FCT a very cost effective convection algorithm. BIC-FCT takes about I5ps 
per point per computational timestep on the CRAY X-MP computer. This is as fast 
as the explicit FCT code currently in use. BIC-FCT has opened up the possibility 
of doing accurate, multidimensional, slow-flow calculations in which fluid expansion 
is important. Because the cost of BIC-FCT is modest even in two dimensions, 
reasonably detailed chemistry models as well as other physical processes can be 
included. Premixed flames, diffusion flames, and turbulent jet flames are some of 
the applications for which BIC-FCT is well suited. 

2.2 Gravity 

Buoyancy effects due to the force of gravity have been incorporated in FLIC. The 
body force term F in Eq. (2.6) is used for this purpose and is given by 

F = 9(P-P~) (2-12) 

where is a suitable reference density, usually the cold ambient density. As 
currently implemented, the direction of the gravity vector is aligned with the flow 
direction, but this restriction can be removed trivially. Indeed, the gravity vector 
can be made time-dependent to simulate g-jitter in microgravity. 

If the ambient density cannot be used, the hydrostatic head has to be included 
in its place. This second approach is not as suitable as the method given by Eq. 
(2.12) because of the need for very high precision calculations. 
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Chemistry 


This chapter describes the integration package TBA that is used to solve the ordi- 
nary differential equations (ODE’s) representing the chemical reactions and energy 
release. TBA is a fully vectorized FORTRAN subroutine for the CRAY X-MP. It is 
designed to replace VSAIM, the older vectorized version of CHEMEQ [8,9], an al- 
gorithm that solves a system of ODE’s in a single computational cell. Both VSAIM 
and TBA are based on CHEMEQ and axe designed to make the calculation of a 
large number of sets of ODE’s more efficient. TBA is faster than VSAIM because 
it is designed to make specific use of the Cray X-MP’s gather /scatter hardware and 
other capabilities. 

The ODE’s solved by these routines are of the form 

^ = F i = Q i -L i n i , (3.1) 

where n,- is the number density of species *, Q, is the rate of formation of species i, 
and T,n, is the rate of destruction of species i. Sometimes these equations can be 
solved by classical algorithms, sometimes they are stiff and need special techniques. 
TBA uses a different algorithm for each type of equation and gains efficiency by 
gathering all equations of a given type together and integrating them by groups. 

A detailed hydrogen-oxygen reaction scheme has been implemented in FLIC. 
This reaction scheme consists of forty-eight reversable reactions involving eight 
species. Nitrogen, acting as a diluent, is considered to be chemically inert. This 
reaction scheme, developed by Burks and Oran [21], has been used by Kailasanath 
et al. [22] in FLAME1D and is given in table 3.1. The toted rate of formation and 
destruction of each species is obtained algebraically from the reaction rates of the 
individual chemical reactions and from the species concentrations. The reaction 
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rate for each chemical reaction r is assumed to follow the modified Arrhenius form: 

k r = AT B e~ c/T (3.2) 

The computation of the rates in Eq. (3.2) is quite time consuming, mainly due 
to the calculation of the exponential term. Some time can be saved by computing 
the temperature dependence of the reaction rates only once per global timestep. 
The rates of formation and destruction of the species, which axe also dependent on 
the species concentrations, is computed as often as needed by TBA, which updates 
them many times during each global timestep. 


3.1 Numerical Method 


TBA uses a second-order predictor-corrector method that is essentially the same 
numerical integration algorithm as CHEMEQ[9]. Normal ODE’s are integrated 
using a simple classical scheme and stiff ODE’s are integrated using an asymptotic 
method. Unlike CHEMEQ or VSAIM, TBA also recognizes equations that will 
approach equilibrium during the period of integration and handles them with a 
third scheme. TBA sorts the equations from every cell into one of three types 
— normal, stiff, or equilibrium, and then integrates all of the equations of each 
type together. A large number of cells are integrated simultaneously; as cells are 
completed, the results are returned to the control program and stored and data 
from new cells are read in and integrated. 

The entire set of equations from each cell is integrated using the smelliest timestep 
required by any equation in the set. The timestep is then increased or decreased 
based on the relative difference between the predictor and corrector stages. The 
predictor stages for the three types of equations are: 


n' = n° x + StFf , 

(Normal) 

(3.3) 

, 0 6tF T 

(Stiff) 

(3.4) 

n' = Qi/Li . 

(Equilibrium) 

(3.5) 

The corrector stages are: 



= < + f [*? + *?] . 

(Normal) 

(3.6) 
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Reaction Rate 

AW 

B 

c< 0) 

Source 

H + HO ^ 0 + H 2 

1.40(— 14) 

1.00 

3.50(4-03) 

[23] 


3.00( — 14) 

1.00 

4.48(4-03) 

[23] 

H + HO 2 ^ H 2 ~h O 2 

4.20( — 11) 

0.00 

3.50(4-02) 

[23] 


9.10( — 11) 

0.00 

2.91(4-04) 

[23] 

H + H0 2 ^ HO + HO 

4.20(— 10) 

0.00 

9.50(4-02) 

[23] 


2.00(— 11) 

0.00 

2.02(4-04) 

[23] 

H 4- H0 2 ^ 0 4- H 2 0 

8.30(— 11) 

0.00 

5.00(4-02) 

[24] 


1.75(— 12) 

0.45 

2.84(4-04) 

k r = Jc,/K c 

H 4- H 2 0 2 ** H0 2 4- H 2 

2.80(— 12) 

0.00 

1.90(4-03) 

[23] 


1.20(— 12) 

0.00 

9.40(4-03) 

[23] 

H 4- H 2 0 2 ** HO + H 2 0 

5.28(— 10) 

0.00 

4.50(4-03) 

[23] 


3.99(— 10) 

0.00 

4.05(4-04) 

kr = k,/K ; 

HO + H 2 ** H 4- H 2 0 

1.83(— 15) 

1.30 

1.84(4-03) 

[25] 


1.79(— 14) 

1.20 

9.61(4-03) 

[25] 

HO + HO ^ H 2 + 0 2 

1.09(— 13) 

0.26 

1.47(4-04) 

II 


2.82(— 11) 

0.00 

2.42(4-04) 

[26] 

HO 4- HO ** 0 + H 2 0 

1.00(— 16) 

1.30 

0.00(4-00) 

[25] 


3.20(— 15) 

1.16 

8.77(4-03) 

kr = kf/K c 

HO 4- H0 2 ^ H 2 0 4- 0 2 

8.30(-ll) 

0.00 

5.03(4-02) 

[27] 


2.38(— 10) 

0.17 

3.69(4-04) 

K = k f /K c 

HO + H 2 0 2 ^ H0 2 4- H 2 

1.70( — 11) 

0.00 

9.10(4-02) 

[23] 


4.70(— 11) 

0.00 

1.65(4-04) 

[23] 

ho 2 + h 2 ^ ho + h 2 o 

1.20( — 12) 

0.00 

9.41(4-03) 

[26] 


1.33(— 14) 

0.43 

3.62(4-04) 

k r = kf/K c 


Table 3.1 Chemical Reaction Rates for H 2 -0 2 Combustion: 

ki = AT B exp( — C/T)^ 

^ Exponentials to the base 10 are given in parentheses: 1.00( — 10) = 1.00 x 10 -10 . 
Bimoleculax reaction rate constants are in units of cm 3 / (molecule s). 
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Reaction Rate 

A («) 

B 

C< a > 

Source 

HO 2 + HO 2 ^ H 2 O 2 4“ O 2 

3.00( — 11) 

0.00 

5.00(4-02) 

[24] 


1.57(-09) 

-0.38 

2.20(4-04) 

k r = kf/K c 

O + HO ^ H + 0 2 

2.72(— 12) 

0.28 

-8.10(4-01) 

kf = k r K c 


3.70(— 10) 

0.00 

8.45(4-03) 

[23] 

0 + H0 2 ^ HO + 0 2 

8.32(— 11) 

0.00 

5.03(4-02) 

(27] 


2.20(— 11) 

0.18 

2.82(4-04) 

k r = kf/K c 

0 + H 2 O 2 ^ H 2 O 4“ O 2 

1.40( — 12) 

0.00 

2.12(4-03) 

[24] 


5.70(— 14) 

0.52 

4.48(4-04) 

K « k f /K c 

O 4- H 2 0 2 ^ HO 4- HO 2 

1.40(— 12) 

0.00 

2.13(4-03) 

[24] 


2.07(— 15) 

0.64 

8.23(4-03) 

kr = kf/K c 

h + h + m^h 2 + m 

1.80(— 30) 

-1.00 

0.00(4-00) 

[23] 


3.70(— 10) 

0.00 

4.83(— 04) 

(23] 

H + HO 4- M ** H 2 0 4- M 

6.20(— 26) 

-2.00 

0.00(4-00) 

[23] 


5.80(— 09) 

0.00 

5.29(4-04) 

[23] 

H + O 2 4- M ^ HO 2 + M 

4.14(— 33) 

0.00 

-5.00(4-02) 

[23] 


3.50(— 09) 

0.00 

2.30(4-04) 

[23] 

HO 4- HO + Mfi H 2 Oa 4- M 

2.50(— 33) 

0.00 

-2.55(4-03) 

[23] 


2.00(— 07) 

0.00 

2.29(4-04) 

[23] 

O + H + M^HO + M 

8.28(-29) 

-1.00 

0.00(4-00) 

[28] 


2.33(— 10) 

0.21 

5.10(4-04) 

k r = kf/Kc 

0 + HO 4- M ^ HO 2 4- M 

2.80(— 31) 

0.00 

0.00(4-00) 

[28] 


1.10(— 04) 

-0.43 

3.22(4-04) 

k r = kf/K, 

o + o + m^o 2 + m 

5.20(— 35) 

0.00 

-9.00(4-02) 

[23] 


3.00(-06) 

-1.00 

5.94(4-04) 

[23] 


Table 3.1 Continued Chemical Reaction Rates for H 2 -O 2 Combustion: k, = 
AT B exp(— C/T)< 6 ) 

^ Exponentials to the base 10 axe given in parentheses: 1.00( — 10) ^ 1.00 x 10 -10 . 
^ Bimolecular reaction rate constants are in units of cm 3 / (molecule s). 
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(Stiff) 


(3.7) 


n. 


n. 


, , g ft (gj ~ jX + jil 
n,+ 4 + ft[L<+L?] 


(Equilibrium) 


(3.8) 


The original CHEMEQ report [9] describes the details of the timestep selection 
algorithm, stiffness criterion, and other important details. A detailed report on 
TBA [29] is currently under preparation. 


3.2 Data- Handling Algorithm 

TBA is designed to handle a large number of cells, each with an independent 
timestep. Thus, each cell is integrated with a different number of subcycles. Cells 
have to be constantly shuffled in and out of the integration routine. Whenever the 
integration of a cell is complete, it is put onto a list. At the end of the integration 
loop, this list is passed to another routine, CDATA, which stores the results and 
inserts new cells to be integrated into the places left by the completed cells. To- 
wards the end of the integration procedure there are no additional cells to integrate 
and the arrays in the integrator will be only partially filled. When this occurs, we 
sort and move all completed cells to the ends of the arrays where they will not be 
integrated further. When the integration of the remaining cells is complete, all the 
data is written out in one operation. Thus we avoid both unnecessary shuffling and 
unnecessary calls to CDATA. The sort algorithm developed specificially for this ap- 
plication is 0(N) and makes the minimum number of swaps. The sort is optimized 
by the use of CRAY assembly routines that make bitwise comparisons. 

In the original VSAIM routine, all equations passed through one integration 
loop where they were tested for stiffness by if statements and either the stiff or 
normal calculations were performed. The CRAY X-MP vectorizes if statements 
by performing the calculations for both possibilities and then throwing away the 
results for the false condition, so this approach is wasteful. TBA creates index 
arrays for different types of equations and sets up a separate integration loop for 
each of the three types. This approach succeeds due to the speed of the CRAY 
X-MP gat her /scat ter hardware. 
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3.3 Programming Strategies 


Many strategies or tricks were used to enhance the speed of TBA, most of which 
are in what is considered extremely poor programming style, but reflect certain 
idiosyncracies of CRAY programming in general. Some are documented in CRAY 
manuals, especially the Optimization Guide [30]. 

For example, many of the arrays are equivalenced as both one- and two-dimen- 
sional arrays. The CRAY will only vectorize the innermost nested loop, so wherever 
possible we looped over the one-dimensional equivalent array, to avoid an outer loop. 
Working space for TBA is supplied in a common block, as this proved substantially 
faster than passing it through the argument list in the calling sequence. Memory 
access is the bottleneck on the CRAY so scalar temporaries are used to reduce 
memory traffic. A full 50% speed up was achieved through their use. A power of 
two as the number of species results in memory conflicts that slow the program 
considerably. There axe eight (2 s ) chemical species in FLIC which would result in 
very poor performance. Thus a “dummy” species was put in, removing the memory 
conflict. 

The following code fragment illustrates the importance of large vector lengths. 
Note that if the order of the loops were switched, the if could be taken out of the 
inner loop and the memory access made contiguous. However, a much shorter inner 
loop results, and the execution time is actually greater. 


do 250 i=l,numeqns 

do 240 j=l,numcells 

if (convchk(j)) then concent(i, j)=corr2d(i, j) 

240 continue 

250 continue 

These sorts of tricks are higly specific to the CRAY X-MP computer and must 
be used to gain full advantage of its speed. Though TBA can be transported to 
other machines (it was developed on a VAX), it was written specifically for the 
CRAY X-MP with its gather/scatter hardware in mind. 
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4. Diffusive Transport Processes 


The diffusive transport processes currently included in FLIC are molecular diffusion 
including the Dufour effect, thermal conduction, and viscosity, ail processes that are 
crucial for describing flame structure. As yet, we have have not included the thermal 
effects on mass diffusion (thermal diffusion or the Soret effect), thermal dissipation 
due to viscosity, or radiation transport. Although it is a second-order effect, thermal 
diffusion may become important for near-limit flames and so should be included in 
the future. Thermal dissipation is negligible at the extremely low Mach numbers 
in the flows under consideration. Radiation effects are important in many flames, 
particularly hydrocarbon flames that form soot, but is not particularly important 
in hydrogen flames. 

The differential equations describing these diffusive processes have been solved 
with a two-dimensional, explicit Eulerian scheme. Spatial derivatives have been ap- 
proximated by central differences and a simple forward- Euler time marching scheme 
has been used for time advancement. This explicit method has a timestep limit 
roughly equal to the timestep required by BIC-FCT for the fluid convection. How- 
ever, each diffusive transport process has its own stability condition, and on oc- 
casion it can require a timestep up to five times smaller. When this is the case, 
the approach we have taken is to subcycle the integration for that process. The 
alternate approach would be to use an implicit algorithm for selected terms, but 
this is generally more expensive than subcycling when fewer than ten subcycles are 
required. 
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4.1 Mass Diffusion 


The part of the species conservation equations for species number density that 
describes mass diffusion are 


= -V ■ (n k V k ) , n k = 1, ...n, p , 

subject to the constraint 

n. p 

E nv; = o . 


(4.1) 


(4.2) 


The effects of molecular diffusion in the energy equation (Dufour effect) appear as 


f = < 4 ' 3 > 

where h k is the temperature-dependent enthalpy for each of the species. The values 
of {Vfc} are found by solving [31] 


Y\V, 

5* = E Tp(V, - H) = VJ? 4 , 

i-1 D *i 


(4.4) 


subject to 

n,p 

E& = o, 

Jt=l 


(4.5) 


and Eq. (4.2), where X k and Y k are the mole and mass fractions of species k. 
The exact solution of this equation for the set {V*} can be obtained by solving the 
matrix equation implied by Eq. ( 4 . 4 ). 


To avoid solving the full matrix equation at each location at each timestep, we 
find an approximation to the {V*} using Fick’s Law and then correct this by the 
procedure described by Coffee and Heimerl [32] to ensure that the constraint Eq. 
(4.2) is met. The diffusion velocity according to Fick’s law is 

R = —Di^VX,, , (4.6) 

where D km is the diffusion coefficient of species k in the mixture of gases. Equation 
4.6 determines the diffusion velocities to within a constant. We then assume that a 
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constant V c is added to all the raw diffusion velocities T4 and require that the sum 
of the diffusion fluxes equal zero. Thus, 


n,p U9p 
£nvt = £n(n+*;) = o 

Jt=l 


n$p 

E 

k=l 


leads to 


K = -£nv i. 

fc— 1 


The component of the corrected diffusion velocity defined by 

V" = v k + v c 


(4.7) 


( 4 . 8 ) 


(4.9) 


is then used in Eq. (4.1). 

The set of {14} found in this way is algebraically equivalent to the first itera- 
tion of the DFLUX algorithm [12], an approach based on a matrix expansion that 
converges to arbitrary order. Values of {14} obtained by the procedure described 
above were compared to the results of DFLUX to check their accuracy. The explicit 
finite differencing used to solve Eq. (4.1) has the numerical stability limit, 

max (D km At/Ax 2 ) < 1/2 . 

Generally the mass-diffusion algorithm is subcycled within an overall timestep de- 
termined by the convection algorithm. However, the code is designed to decrease 
the overall timestep to below that required by the mass-diffusion stability limit if 
the number of subcycles required exceeds a specified maximum value. Subcycling 
becomes necessary when the temperature of the reacting flow becomes high and the 
diffusion coefficients increase accordingly. 

Determining the diffusion velocities by Eq. (4.6) requires as input the set of dif- 
fusion coefficients of species k into the mixture, {Z?* m }. However, these quantities 
are difficult to obtain from first principles and axe usually found by applying a mix- 
ture rule to the individual binary diffusion coefficients. Binary diffusion coefficients 
can be estimated theoretically [33] and sometimes measured experimentally [33,34]. 
Here we use the same approach as in FLAMElD (Kailasanath et al. [22]). Binary 
diffusion coefficients are expressed in the form 

D m = — T Bk ‘ ( 4 . 10 ) 

ntot 
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0 

h 2 

OH 

h 2 o 

02 

H0 2 

H 2 0 2 

n 2 

H 

6.30(17) 

8.29(17) 

6.30(17) 

6.70(17) 

6.70(17) 

6.70(17) 

4.43(17) 

6.10(17) 


7.28(-l) 

7.28(-l) 

7.28(-l) 

7.28(-l) 

7.32(-l) 

7.32(-l) 

7.28(-l) 

7.32(-l) 

0 


3.61(17) 

1.22(17) 

2.73(17) 

9.69(16) 

9.69(16) 

1.57(17) 

9.69(16) 



7.32(-l) 

7.74(-l) 

6.32(-l) 

7.74(-l) 

7.74(-l) 

6.32(-l) 

7.74(-l) 

h 2 



3.49(17) 

6.41(17) 

3.06(17) 

3.06(17) 

4.02(17) 

2.84(17) 




7.32(-l) 

6.32(-l) 

7.32(-l) 

7.32(-l) 

6.32(-l) 

7.38(-l) 

OH 




2.73(17) 

1.16(17) 

9.69(16) 

1.57(17) 

9.69(16) 





6.32(-l) 

7.24(-l) 

7.74(-l) 

6.32(-l) 

7.74(-l) 

h 2 o 





2.04(17) 

2.04(17) 

1.57(17) 

1.89(17) 






6.32(-l) 

6.32(-l) 

6.32(-l) 

6.32(-l) 

o 2 






8.74(17) 

1.14(17) 

8.29(16) 







7.24(-l) 

6.32(-l) 

7.24(-l) 

ho 2 







1.14(17) 

8.85(16) 








6.32(-l) 

7.74(-l) 

h 2 o 2 








1.14(17) 









6.32(-l) 


T B i k 

Table 4.1 Binary diffusion coefficients expressed in the form: Dj k = A ]k — . For 
each pair of species, the upper term is Aj k and the lower term is B jk , exponentials 
to the base 10 are given in parentheses. 


where the sets {A k i} and {B k i} are tabulated [22] for each pair of species in the 
hydrogen-oxygen reaction system, and are given in table 4.1. These are then com- 
bined by a mixture rule [35,36] 


Dkm 


i-n 

^ Xi_ ’ 

t(Du 


(4.11) 


which provides values of D km to use in Eq. (4.6). 
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4.2 Thermal Conduction 


The effects of thermal conduction axe expressed in the energy equation as 


dE 

dt 


-V . (kVT) , 


(4.12) 


where k is the mixture thermal conductivity. Explicit finite differencing introduces 
a stability limit 

max (nAt/pCyAx 2 } < 1/2 , 

where n/pCp is the thermal diffusivity coefficient. Subcycling is used here in the 
same manner as it is for mass diffusion. However, this stability condition is less 
stringent than that for mass diffusion and typically only two or three subcycles are 
needed. 


The mixture thermal conductivity k is obtained by combining the thermal con- 
ductivities of the individual gases {«*} that are in the mixture. The {«*} are 
estimated theoretically and are a function of temperature. We have used the third- 
order polynomial fit in temperature determined by Laskey [37] and is presented in 
Table 4.2. The mixture thermal conductivity is then calculated using (Mathur et 


al. [38]) 


K 




* 

1 

Tlgp 

E + 

fc= 1 

1 

2 




fc=l K k . 


(4.13) 


4.3 Viscosity 


The viscosity terms in the Navier-Stokes equations are included in FLIC. The 
stress- tensor term, which represents the effect of viscosity in the momentum con- 
servation equations, is: 


dpV 

dt 


-V-r , 


(4.14) 


where 

I 


(4.15) 
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Species 

A B C D 

H 

4.710(3) 3.354(1) -9.971(-3) 1.964(-6) 

0 

1.089(3) 1.038(1) -3.739(-3) 8.251(-7) 

h 2 

6.306(3) 4.304(1) -8.505(-3) 2.160(-6) 

OH 

1.679(3) 1.091(1) -1.613(-3) 4.150(-7) 

H 2 0 

-2.077(2) 1.603(1) -7.932(-4) 1.530(-7) 

o 2 

3.862(2) 8.613(0) -1.966(-3) 3.619(-7) 

ho 2 

1.576(2) 1.070(1) -1.143(-3) 4.471(-8) 

h 2 o 2 

-1.097(3) 9.895(0) -1.779(-3) 1.396(-7) 

n 2 

7.024(2) 6.917(0) -1.191(-3) 2.035(-7) 


Table 4.2 Thermal conductivity of species k, = A + BT + CT 2 + DT 3 , 
erg/cm-s-K. Exponentials to the base 10 are given in parentheses. 

and p is the dynamic viscosity coefficient. The quantity A, the second coefficient of 
viscosity, is set to —2/Zp. The stress tensor r includes all the viscous terms which 
arise in the compressible Navier-Stokes equations. 

Eq. (4.14) is solved explicitly in the same manner as the mass diffusion or 
thermal conduction equations. Thus there is a stability criterion given by 

max (pAt/pAx 2 ') < 1/2 , 

where p/p is the kinematic viscosity. Subcycling is used so that the overall com- 
putational timestep is set by the convection stability limit and not by the viscosity 
stability limit. If the number of subcycles required for stability exceeds some max- 
imum value, the global timestep is reduced. The viscous diffusion algorithm was 
tested using two test problems: 1) the boundary-layer growth over a flat plate par- 
allel to the flow, and 2) the boundary-layer thickness on a flat plate normal to the 
flow (stagnation point flow). 

For the parallel-plate test, the velocity profile of the parallel flow was initially 
uniform along the plate with a typical velocity profile that is valid for a boundary- 
layer thickness greater than three computational cells. This particular initialization 
was chosen because a physical boundary layer less than three cells wide would be 
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swamped by numerical diffusion from the FCT algorithm that creates a boundary 
layer at least this thick. Setting up the problem this way simulates the growth 
of a boundary layer away from the leading edge of the plate and minimizes any 
numerical effects on the solution. The result of this test is a boundary layer whose 
growth matches the Blasius solution. 

The stagnation-flow test was initialized with a boundary layer of constant thick- 
ness and uniform velocity along the plate. Again, the initial boundary layer thick- 
ness was more than three cells to minimize numerical effects. The results showed 
the development of a constant-thickness boundary layer whose velocity profile very 
closely matched that predicted by theory. 


For a gas containing a single species k, the dynamic viscosity fik can be otained 
from the kinetic theory of gases [33]. Over a suitable range of temperature, this can 
be expressed as a third-order polynomial in temperature. Laskey [37] has compared 
this polynomial fit, presented in Table 4.3, to the calculations and to tabulated 
values for the viscosity and found good agreement. The mixture dynamic viscosity 
is calculated using the expression (Wilke [39]) 


/* — 2J —p > ( 4 - 16 ) 

* =1 £*;**i 


where 



(4.17) 


This mixture viscosity, /i, is used in the stress tensor (Eq. (4.15)) to model the 
viscous portion of the Navier-Stokes equations. 
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Species 

A 

B 

C 

D 

H 

1.516(-5) 

1.074(-7) 

-3.178(-11) 

6.255(-15) 

0 

4.504(-5) 

5.355(-7) 

-1.811(-10) 

3.740(-14) 

h 2 

2.802(-5) 

2.236(-7) 

-6.958(-ll) 

1.399(-14) 

OH 

5.630(-5) 

5.193(-7) 

-1.678(-10) 

3.413(-14) 

H 2 0 

-8.597(-6) 

6.608(-7) 

-2.305(-10) 

4.601(-14) 

o 2 

4.930(-5) 

5.861(-7) 

-1.983(-10) 

4.093(-14) 

ho 2 

5.006(-5) 

5.952(-7) 

-2.013(-10) 

4.156(-14) 

h 2 o 2 

-9.109(-6) 

3.966(-7) 

-1.345(-10) 

2.623(-14) 

n 2 

5.302(-5) 

4.596(-7) 

-1.464(-10) 

2.969(-14) 


Table 4.3 Viscosity of species k,fi k = A + BT + CT 2 + DT 3 , dyne-s/cm 2 . Expo- 
nentials to the base 10 are given in parentheses. 
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5. Model Integration 


The conservation equations contain terms representing convection, buoyancy, ther- 
mal conduction, molecular diffusion, viscous diffusion, and chemical reactions. The 
approach that FLIC uses is to determine a global timestep, solve equations rep- 
resenting the individual physical processes separately for that timestep, and then 
couple the solutions. The coupling procedure is a variation of the standard timestep 
splitting method of Yanenko[40] and described for reactive flows by Oran and 
Boris [15]Chapters 4 and 13. 

The usual explicit timestep-splitting approach assumes that in some predeter- 
mined global timestep, the effect of all the physical processes can be evaluated as 
a running sum of the effects of individual processes. Each physical process is inte- 
grated independently using the results of the previous process as initial conditions. 
This method is correct in the limit of small timesteps and works well in a practical 
sense when the changes in the variables during the global timestep axe small. Using 
this appproach, the global timestep is often limited to the smallest timestep required 
by the stability limits of the integration algorithms for the various processes. This 
is the approach we have used in a number of programs in which the convection is 
solved by an explicit integration procedure. The global timestep is usually deter- 
mined by the CFL condition on the sum of the sound speed and the fluid velocity. 
The chemistry integration is subcycled in the global timestep. However, subcycling 
can sometimes be used for individual processes if changes in variables due to that 
particular process are not too large. 

Figure 5.1 summarizes the integration process in a typical FLIC timestep. The 
global timestep is first estimated based on the stability requirements of the con- 
vection algorithm, BIC-FCT. Then it is compared with the stability requirements 
of the particular physical processes which are allowed to subcycle if necesary to 
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Figure 5,1 Flowchart of FLIC, indicating interprocess communication. 
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ensure stability. The overall timestep must be decreased sometimes to ensure that 
a physical variable does not change too much during a timestep. The most obvious 
difference between the scheme shown in Fig. 5.1 and the standard timestep splitting 
approach is that changes in the internal energy caused by each process are evaluated 
and added up at the end of a global timestep. This energy change is then used by 
the BIC-FCT convection algorithm (see Eq. (2.7)). This approach is quite similar 
to that used in FLAME1D [22], however, now changes in the internal energy are 
accumulated instead of changes in the pressure. The ADINC method [11], used in 
FLAME1D, does not solve an energy equation, and thus the only way other pro- 
cesses cam be coupled is through the pressure. The BIC-FCT scheme used by FLIC 
does not require that energy changes be applied only during the convection step; 
this is merely a convenience that allows for larger global timesteps due to tighter 
coupling. 

Some of the individual process integrations in FLIC are subcycled within a 
global timestep, including the ordinary differential equations representing the chem- 
ical reactions and the diffusive terms such as molecular diffusion and thermal con- 
duction. For example, the timestep limit imposed by some chemical reactions may 
be orders of magnitude lower than that required by other physical processes, and 
so the chemistry integration is subcycled. Subcycling is built into the chemistry in- 
tegration in an extremely sophisticated manner [9], so that the maximum allowable 
timestep at each computational cell is used, completely independent of the timestep 
in other cells. The chemistry is integrated up to the overall timestep before it is 
coupled to the other processes. However, if the energy release in an overall timestep 
due to chemistry is too large (typically greater than 10%), then the overall timestep 
must be decreased. Mass diffusion and thermal conduction are also subcycled, but 
only up to five times. The accuracy of the solution is generally tested by performing 
a separate calculation with a smaller timestep and noting whether the solution has 
converged. 

All dependent variables, except for internal energy, are updated after each pro- 
cess integration, but dependent variables are not updated during subcycling of a 
process. This leads to a considerable savings, especially during the chemistry in- 
tegration, because the evaluation of the temperature exponentials in Arrhenius ex- 
pressions is done only once per global timestep. On the other hand, the global 
timestep may have to be decreased if the dependent variables change too much. 
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One consequence of a source term in the energy equation is the possibility of 
“ripples” in the pressure, which then manifest themselves in other variables as well. 
These ripples arise if a strong source is localized to a region only a few cells wide. 
These ripples are insignificant in FLIC where the flame zone is well resolved, but 
can be significant in other applications [37]. One way to avoid the ripples is to use 
a high-frequency filter such as 

p/i/tererf = p + oV 4j> 


where a is a small constant. Other filters, including another FCT step, can also be 
used. 

The particular advantages to using timestep splitting are that we can write very 
modular programs in which the integration of each physical process can be carried 
out with an optimum method, debugging is simpler, and explicit subcycling can 
be used to keep the costs down. Implicit methods can be used to avoid subcy- 
cling, but are more expensive when compared to explicit methods subcycled only 
a few times. In FLIC, only five subcycles are allowed for each of the diffusion 
processes. Chemistry, on the other hand, subcycles thousands of times. At this 
point, it is not entirely clear if the extreme simplicity of the CHEMEQ scheme [9] 
results in faster integration of the chemistry equations than a more complicated 
implicit scheme. Disadvantages of timestep splitting are that the coupling process 
can be complicated, the algorithms and the overall program are less stable, and the 
timestep must be carefully controlled. However, we have found that the benefits of 
modular and fast programs outweigh potential disadvantages. This is discussed in 
some detail by Oran and Boris [15], pages 131-133. 
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6. Applications 


FLIC has been written in a general manner so that it can be applied to solve 
a variety of slowly evolving problems involving chemistry and diffusive transport 
processes. It has already been applied to the study of multidimensional flames in 
premixed gases [1,41,42] and co-flowing diffusion flames [37,43]. A specific applica- 
tion of the code to the study of cellular flames in hydrogen-oxygen-nitrogen mixtures 
is described below to show what is required to perform a calculation with the code 
and to interpret the output from the code. Then other applications of the code jure 
briefly discussed to bring out the generality of FLIC. 

6.1 A Sample Calculation 

Flames in lean hydrogen-oxygen-nitrogen mixtures are known to exhibit a multidi- 
mensional cellular structure. The cellular structure is the result of a thermo-diffusive 
instability of a planar flame in the same mixture. Below we demonstrate how FLIC 
can be used to simulate the transition from a planar flame to a multidimensional 
cellular flame in a zero-gravity environment. 

We need as input to the model: 

A chemical reaction scheme involving all the species of interest (table 3.1), 
Molecular diffusion coefficients for each pair of species (table 4.1), 

For each species: 

Molecular weights, 

Thermal conductivity (Table 4.2), 

Viscosity (Table 4.3), 
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Heats of formation [22], 
Enthalpy coefficients [22]. 


To complete the specification of the problem, we also need the initial and bound- 
ary conditions such as composition, pressure, and temperature, in addition to the 
physical and chemical parameters given above. Figure 6.1 describes the configu- 
ration studied and gives the boundary conditions of the computational domain. 
Unburnt gas flows in from the left, and reaction products of the flame front flow out 
at the right. If the inlet velocity equals the burning velocity of the flame, the flame 
is fixed in space yielding a steady solution. Thus, transient effects from ignition 
can be eliminated. The initial conditions for the two-dimensional calculations were 
taken from one-dimensional calculations that gave the conditions for steady, propa- 
gating flames. The two-dimensional computational domain for this simulation was 
2.0 cm x 4.5 cm, resolved by a 56 x 96 variably spaced grid. Fine zones, 0.36 mm 
x 0.15 mm, were clustered around the flame front. 

The initial conditions specify a planar flame in a fuel-lean hydrogen-oxygen 
mixture diluted with nitrogen, HjiOjiNj/l.SrlrlO, aflame that showed multidimen- 
sional structure in the experiments by Mitani and Williams [4]. In order to study 
the evolution to cellular structure, the initial conditions were perturbed by displac- 
ing the center portion of the planar flame in the direction of the flow. The evolution 
to cellular structure is obtained by studying the output from the simulations. 

The output from the calculation consists of the spatial and temporal distribution 
of all the species involved as well as the fluid density, temperature, pressure, internal 
energy, and momenta. Display of the data is not done by FLIC , but is instead done 
as a post-processing operation. This cuts down the length and complexity of the 
FLIC code. A very useful diagnostic is contour plots showing isotherms and species 
distributions. Diagnostics, such as color-flood plots, allow visual interpretation of 
the data. 

Figure 6.2 shows a sequence of isotherms from the calculation. The isotherms 
show that the temperature increases in the center portion of the flame, convex to the 
flow, and decreases in the the two adjacent concave regions, indicating more vigorous 
reaction in the convex region. The atomic hydrogen concentration increases in the 
convex and decreases in the concave regions. Also, the burning velocity in the 
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Figure 6.1 Initial and boundary conditions for the two-dimensional flame calcu- 
lations. 


30 


A-34 




convex region appears slightly higher than the burning velocity in a planar region, 
and the burning velocity in the concave regions were noticeably lower. Thus, this 
calculation shows that the planar lean one-dimensional hydrogen flame is unstable 
and evolves into a multidimensional flame having a cellular structure. 


6.2 Applications of FLIC 


The FLIC code has been used extensively to study the detailed evolution of cel- 
lular flames in premixed gases [1], to investigate the mechanisms which can lead to 
cellular structure [42], effects of gravity on flame structure and instabilities [41,42]. 
These studies have helped strengthen the prevailing theory of cellular instability 
and cast serious doubt on another theory which was also under consideration [1] 
FLIC has been geared primarily to these types of applications and several other 
related applications to premixed flames are planned for the near future. 

FLIC can be converted to the study of diffusion flames extremely easily, and 
is begging for a suitable application in this area. A low-speed diffusion flame 
code [37,43] which has been used to study jet flames uses the same transport and 
diffusive packages as FLIC. This code uses simplified chemistry, which was first 
calibrated by comparison to a detailed calculation with FLIC. 


6.3 Future Applications and Code Development 


Calculations such as the one described earlier are time consuming (5 hours of 
CRAY X-MP time). There is interest to carry out these calculations in a larger 
domain for longer times. The cellular structure exhibited by flames is actually three- 
dimensional, so if these instabilities are to be studied fully, a three-dimensional ver- 
sion of FLIC is required. The computer time for these studies can quickly become 
intolerable; it is proposed to alleviate this by taking advantage of the multiprocessor 
architechure of the CRAY X-MP and the new Y-MP with its 16 processors. This 
will require some restructuring of the code and its algorithms to utilize microtask- 
ing appropriately. This restructuring may be simplified with the new “autotasking” 
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facility on the CRAY Y-MP. 


One area of interest is the investigation of the behavior of flames near extinguish- 
ment. In particular, the prediction of flammability limits of flames in a flammability 
tube will provide a means for quantitative comparison with experiment. This com- 
plex task will require the addition of mechanisms which represent losses, including 
losses to the wail. The physics of the loss mechanisms of radicals to the walls is not 
yet fully understood. These additional mechanisms will need to be incorporated 
into FLIC. 

Detailed calculations will need to be performed for other more complex fuels 
of practical interest. The primary difficulty is in coping with the large number of 
species and chemical reactions that will be required for even the simplest hydrocar- 
bon fuels. While the precise scaling of computer time with chemical complexity is 
not known, it is expected that the increase in the number of species will have the 
most dramatic effect. Thus, there is a need for reliable simplified chemistry models 
which have been first validated against a full model. It is anticipated that suitable 
simplified models for methane will become available soon. Radiation can not be 
neglected in hydrocarbon combustion. Thus a gas phase radiation model will have 
to be incorporated. Consideration of soot will have to wait until reliable models are 
available, which will not be in the near future. 

All future applications are in areas that will require enormous amounts of com- 
puter time. One way to cut down on costs is to perform calculations only where 
they axe strictly required. Calculations can be skipped in regions of low gradients, 
be it spatial or temporal gradients. This will require dynamic regridding or, more 
generally, re-discretization of the governing equations. A suitably efficient algorithm 
for this is essential, and its development is underway. Limitations and peculiarities 
of the target computer architecture will play a large role in determining the shape 
of these future versions of FLIC. 
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Abstract 

A time-dependent, one-dlaenmlonml, Lagranglan model 
was used to study laminar f laaoo la stoicbloaotrlc hydrogen- 
oxygaa mixtures dllatad with altrogaa. For atolcbloastrlc 
bydrogaa-alr mixtures wo have aaaa that a spherically ex- 
panding flaae first dacolsratos oatil tho velocity reaches 
a alnlma value, aad then it accelerates. For large radii, 
the burning velocity approaches tho planar burning veloc- 
ity. These saao trends are also observed as the aaount 
of diluent is Increased. With increasing dilution, the 
flaaes reach their alnlma velocities at later tlaes and 
larger radii. These observations are explained on the 
basis of flaao stretch. The spherical geoaetry results 
are coapared with the results froa another set of calcu- 
lations in planar geoaetry. These show that the alniauu 
burning velocity reached by a spherical flaae is less than 
that of a planar flaae in the saae alxture. In both pla- 
nar and spherical geoaetrles. the effect of Increasing the 
dilution is to lower the burning velocities. Since the 
burning velocity is smaller in the spherical geometry, the 
flaae can be extinguished (or quenched) with less dilution 
in the spherical geometry than in the planar geometry. We 
also discuss the implications of these results to laminar 
flaae quenching and flammability limits. 
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Introduction 

In this paper, we present and discuss time -dependent 
calculations of one-dimensional laminar flames in stoi- 
chiometric hydrogen- oxygen mixtures diluted with nitrogen. 
We use these calculations to study the effects of cur- 
vature Cor stretch) and dilution on flame propagation in 
preaixed gases. The results presented also have important 
applications to laminar flame quenching and flammability 
limits in the absence of external heat sinks. 

Previous work on propane-air sixtures by Strehlow 1 
shows that rich mixtures exhibit a higher flame velocity 
at small radii than at large radii, while lean mixtures 
show a lower flame velocity at small radii than at large 
radii. In all cases, for large enough radii the flame ve- 
locity relaxes to the appropriate planar flame velocity. 
Strehlow attributes these effects to preferential diffu- 
sion of the lighter species, and this has been verified 
by Frankel and Shlvaahlnsky 2 using an asymptotic analy- 
sis. More recently. Law 3 reused an asymptotic analysis to 
explain the behavior of stretched flames in rich and lean 
mixtures of nethane and propane is air. These theorotical 
analyses are for either fuel-rich or fuel-lean conditions 
and use a simple phenomenology to represent the chemical 
kinetics. 

This paper presents calculations of planar and spheri- 
cally expanding flanes in stoichiometric hydrogen-oxygen 
mixtures diluted with nitrogen. The results presented 
here are obtained from numerical simulations, which use a 
time-dependent, one -dimensional. Lagrangian model s . This 
numerical model was developed specifically to study the 
initiation, propagation and quenching of laminar flames. 

In addition to using a detailed chemical reaction mecha- 
nism, the model includes the effects of molecular diffu- 
sion, thermal conduction, and thermal diffusion of the in- 
dividual species considered. Because of the level of de- 
tail incorporated in the model, the spatial structure and 
temporal evolution of the flame structure can be highly 
resolved. 


The Numerical Model 

The numerical model solves the time-dependent conser- 
vation equations for mass, momentum, and energy 6 > 7 .- The 
model has been used for a variety of flame studies, in- 
cluding calculations of minimum ignition energies 8 , 9 quench 
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volumes 8 , and burning velocities 1 °: The model has a mod- 

ular form and persits a wide variety of geometric, ini- 
tial. boundary, and time varying energy input conditions. 
The algorithms representing the various chemical and phys- 
ical processes are integrated separately and then asymp- 
totically coupled by time-step splitting techniques 7 . The 
convective transport is solved by the algorithm ADINC, a 
Lagrangian fluid dynamic algorithm that solves implicitly 
for the pressures 1 K The method gives an accurate repre- 
sentation of material interfaces and allows steep gradi- 
ents in species and temperature to develop and be main- 
tained. In addition to considering the thermal conduction 
and molecular diffusion processes in detail, the model 
also includes thermal diffusion. The chemical interac- 
tions are described by a set of nonlinear, coupled ordi- 
nary differential equations that are solved using a fully 
vectorized version of the selected asymptotic integration 
method CHEMEQ 1 2 A *. 

For this study, we have used the hydrogen- oxygen re- 
action scheme 14 , which Involves the eight reactive species 
H> , 0,. H. 0. QR, HjO. HO), R a 0) and a diluent, chosen to 
be nitrogen. The thermochenlcal properties of the var- 
ious species involved are taken from the JAHAF tables 1 **- 
The chemical reaction mechanism (given in Kailasanath et 
al. V ) has been extensively tested and shown to give good 
results. Burks and Oran 1 5 showed that the results com- 
puted with this mechanism compare well with experimentally 
observed induction times, second explosion limits, and the 

i 



Flg.l The flow velocity end teaperature profiles In e planar 
fleae propagating in a Hj:0i:Nj/2:1:4 mixture. 
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Fig. 2 The flow velocity ud temperature prof 11* a at a partic- 
ular tlaa in a spherically axpanding fleas iu a Hj : Qj :Hi/2: 1 : 4 alxtura. 


tnporil behavior of reactive species. Oran at al. i6 have 
•ton that tha aachanlan coupled to • convective trans- 
port algorithm gives good results in the slsulatlon of the 
conditions behind s reflected shock. The laalnar burning 
velocities calculated using the aechanisa are in agreement 
with experinental data 1 °* 

( For the calculations presented below, the model was 
configured with an open boundary at one end to simulate an 
unconfined system. Most of the calculations were done in 
a spherically symmetric one -dimensional geometry. Other 
calculations were done assuming a planar configuration. 

All the calculations are for stoichiometric mixtures of 
hydrogen and oxygen at an initial temperature and pressure 
of 208 K and 1 atm, respectively. The amount of nitrogen 
was varied, thereby varying the dilution. 

Results and Discussion 
Estimation of Burning Velocity 

For either thin or planar flames, the instantaneous 
normal burning velocity can be calculated from the flame 
velocity if we know the velocity of the unbumt gases ahead 
of the flame. For planar flames, the velocity of the un- 
burnt gases ahead of the flame is constant, as shown in 
Fig. 1, where the spatial variation of the flow veloc- 
ity and the temperature across the flame is shown for a 
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hydrogen-air mixture. Hence the burning velocity can be 
unambiguously determined as the difference between the 
flame velocity and the flow velocity ahead of the flame, 

V b = V f - V flow 


For a thin flame, a similar definition ia adequate. How- 
ever, when the flame has a finite thickness and is curved 
(as for a spherically expanding flame) , the appropriate 
definitions for the location of the flame front and the 
fluid velocity of the unburnt gases are ambiguous. This 
can be seen in Fig. 2, in which the spatial variation of 
the flow velocity and the temperature across a spheri- 
cal flams in a hydrogen-air mixture is shown. The fluid 
velocity reaches a maximum within the flame and then de- 
creases ahead of the flame. 

For the flames studied in. this paper, two reference 
fluid velocities have been chosen for estimating the burn- 
ing velocity. One of them is the maximum fluid velocity 
Vmam in the system, and the other is the fluid velocity of 
the unburnt gases V 300 corresponding to the first location 
ahead of the flams with temperature of 300 K. The lower 
estimate for the burning velocity Is obtained as V/ - V max 
and the upper estimate as V/ - V SO o . For a planar flame 
in which the fluid velocity ahead of the flame is constant 
and the sane as the maximum fluid velocity, the two esti- 
mates for the burning velocity are identical. 

1 

Effects of Curvature 

In the first set of calculations, a spherically ex- 
panding flame in a stoichiometric hydrogen-air (actually, 
H3:02:N)/2:1:4) mixture was studied. In these simula- 
tions, energy was deposited linearly over a fixed period 
of time at the center of a spherically symmetric system. 

The radius of energy deposition was larger than the quench- 
radlus 9 and was held constant. We then tracked the spa- 
tial location of the flame kernel as a function of time 
and used this to calculate the apparent flame velocity 
as a function of time. The results of such a calculation 
are shown in Fig. 3. The flame velocity Vj Initially de- 
creases with time until it reaches a minimum value, and 
then it Increases. The figure also shows the maximum fluid 
velocity V m am of the system and the fluid velocity of the 
unburnt gases Vjoo corresponding to the first location 
ahead of the flame with temperature of 300 K. With in- 
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creasing radii, the difference between the two fluid ve- 
locities decreases. For very large radii (not shown in 
the figure), there is an unambiguous burning velocity, 
since the two fluid velocities approach the same value. 
This estimated burning velocity approaches both the ex- 
perimental value and the value determined from a separate 
planar one-dimensional calculation. 

SXXtcf oX flllafclgn 

We also performed a series of calculations in which 
the amount of diluent was varied. The results of one such 
calculation, for the H):0]:Nj/2:l:7 mixture, are shown in 
Fig. 4. Increasing the amount of diluent does not change 
the observed trend*: The flame velocity and the burning 

velocity decrease with Increasing radii, attain a mini- 
mum value, and then increase with increasing radii. This 
calculation and the results shown in Fig. 5 for the 2:1:10 
case show that it takes longer for the flame to reach the 
minimum flame velocity as dilution increases. For the 
2:1:10 mixture, we expect the flame velocity to increase 
again at larger radius, following the trend seen in the 
other two mixtures. For this mixture, the burning veloc- 
ity corresponding to the minimum flame velocity is 0.30 
- 0.35 s/s, whereas the calculated burning velocity for a 



Fig. 3 Time history of the prop- 
agation of a spherically expanding 
flame in a H 2 :02 : N 2 / 2 : 1 :4 mixture. 
The flame velocity is denoted Vj% 
the maximum fluid velocity is 
V maa , and the velocity of the 
first position ahead of the flame 
with a temperature of 300 K is 
^soo- 
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Fig. 4 Time history of the 
propagation of a spherically 
expanding flame in H 2 : 02 :N 2 / 
2:1:7 mixture. (See Fig. 3 
for legends.) 


planar liana in the aaaa mixture is 0.85 - 0.00 a/s. 

Thus tha burning velocity of a spherical liana correspond- 
ing to tha ainiaun valua ol tha llaaa velocity is smaller 
than that of a planar llaaa in tha saaa mixture. We have 
observed this trend in all ol the stoichiometric mixtures 
we have studied. 

Soaa ol tha more dilute mixtures studied do not sup- 
port spherical llaaa propagation. One such case is the 
Hj ^0j :N a /2: 1 : 13 mixture, in which tha flame velocity does 
not level oil but continues to decrease until the flame 
dies. Figure 0 shows V/, V maa . and V 30 o as a function ol 
tine lor this case. We varied tha radius ol energy depo- 
sition, tha amount ol energy input, and the mode ol energy 
input. In all cases, the llaaa died alter propagating for 
a short time. 

When wa modeled llama propagation in the planar geom- 
etry, we observed signillciant burning velocities lor a 
wider range ol mixtures than in the spherical geometry. 

For example, the burning velocity was between 0.34 and 
0.40 a/s lor the H3:0):N)/2:1:13 mixture. As mentioned 
above, this mixture did not support spherical flame propa- 
gation. 

We have used the term "spherical flame" to mean a flame 
whose radius ol curvature is comparable to its thickness. 

A llama of very large radius (typically, an order of mag- 
nitude larger than the flame thickness) behaves like a 
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TIME W 


Fig. 6 Time history oi tho propagation of • spherically expanding 
flase In Ht:0t:Vj/3: 1: 10 mixture. (So# Fig. 3 for logonds.) 


planar liana, and can propagate in the Hj : Oj : Nj/2: 1 : 13 
mixture . Such a liana can bn initiated by dopoaiting a 
larga amount ol energy over a vary larga spherical volume 
oi gas. 

Effects ol Stretch 

Tha observed tine histories of the liana velocity can 
be explained on the basis ol stretch (due to flow diver- 
gence) and eheaical kinetics. First, consider the initial 
deceleration ol V/ in spherical geonetry. Initially, due 
to How divergence, the energy released in chenlcal reac- 
tions does not balance the energy conducted and dillused 
into the unreacted mixture. Because ol this, the flame 
velocity and temperature decrease as the flame expands. 

The decrease in the liana temperature decreases the energy 
release rate, which leads to a lower llame velocity, and 
so on. However, this process does not continue indefi- 
nitely. because the stretch ellect decreases with increas- 
ing radii. Because ol this ellect a minimum llame veloc- 
ity is observed in spherical llame propagation. As the 
radius ol the llama Increases lurther, the energy released 
in chemical reactions is larger than the en- 
ergy conducted and dillused into the unburnt mixture. Now 
the llame velocity Increases with increasing radii until 
a balance is attained between the energy released by chem- 
ical reactions and the energy conducted and diffused into 
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TIME (ms> 

Fif.6 Ti** history of tho propagation of t spherically expanding 
flu* la a Hi :0j :H)/2: 1 : 13 mixture. (Soa Fig. 3 for legends.) 



120 160 200 240 280 320 

TIME liitl 

Tig. 7 Tlaa history of tho aaxlana temperature of a spherically 
expanding flaaa In a H t :0,:I,/2:l:4 alxtnra. 


the unburnt alxtnra ahead of tha Haas. Such a balance 
occur* lor large radii, and than tha flaae propagates as 
if it vara planar. 

Tha above observations are verified by the tempera- 
tures calculated in the numerical simulations. The time 
history of tha aaxlnua temperature in a spherically ex- 
panding flame in a hydrogen-air mixture is shown in Fig. 7. 
Tha maximum temperature corresponding to the minimum ve- 
locity spherical flame is significantly lower than that of 
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a planar flame in the same mixture. After reaching a min- 
imum, the flame temperature increases and tends toward the 
temperature corresponding to a planar flame in the same 
mixture . 

The actual radius at which the minimum flame velocity 
occurs and the radius at which the effect of stretch be- 
comes negligible depend on the flame thicknesss. Since 
the flame thickness increases with increasing dilution, 
a flame in a more dilute mixture will have to propagate 
to a larger radius in order to overcome the effects of 
stretch. This explains the observation made above, which 
showed that the minimum flame velocity is reached at a 
later time (and a larger radius) as dilution is increased. 
Not only is the minimum flame velocity reached later, but 
the magnitude of the velocity and the flame temperature 
are smaller. Therefore, by increasing the dilution, one 
can obtain a mixture which will be quenched by endothermic 
reactions before it can overcome the effects of stretch. 
This appears to be the case with the H 2 :0 a :N 2 /2:l:13 mix- 
ture. 


Summary and Conclusions 

We have studied flames in stoichiometric hydrogen- 
oxygen mixtures diluted by nitrogen in spherical and pla- 
nar geometries using time -dependent , one -dimensional nu- 
merical simulations. The model includes a full descrip- 
tion of the eight reactive species involved in hydrogen- 
oxygen kiitetics, plus detailed models for the chemical ki- 
netics, thermal conduction, and thermal and molecular dif- 
fusion processes. These studies have led to a number of 
interesting observations and important questions on flame 
behavior. 

We have seen that in a stoichiometric hydrogen-air 
mixture, the flame velocity of spherically expanding flames 
first decreases with increasing radii and then increases. 
For large radii, the burning velocity approaches the pla- 
nar burning velocity. These same trends are observed when 
the amount of diluent is increased. However, with in- 
crease in dilution, the minimus flame velocity is reached 
at a later time and at a larger radius. Comparing these 
spherical geometry results with another set of calcula- 
tions in planar geometry shows that the minimum burning 
velocity of a spherical flame is less than that of a pla- 
nar flame in the same mixture. 
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In both planar and spherical geometries, the effect of 
increasing the dilution is to lower the burning veloci- 
ties. Since the burning velocity is smaller in the spher- 
ical geometry, the flame can be extinguished (or quenched) 
with less dilution in the spherical geometry than in the 
planar geometry. The existence of a quench-radius and a 
minimum ignition energy in the absence of heat or radi- 
cal loss to walls or other external sinks has already been 
demonstrated 9 . The results presented here indicate that 
such "self -quenching" also depends on the geometry or cur- 
vature . 

The time histories of the flame velocity are controlled 
by stretch (due to flow divergence) and energy release by 
chemical reactions. The maximum temperature of a spher- 
ically expanding flame exhibits the same trends as the 
flame velocity. The actual radius at which the minimum 
flame velocity occurs and the radius at which the effect 
of stretch becomes less important depends on the thickness 
of the flams. The thickness ol the flame increases with 
increasing dilution. Therefore, a more dilute flame will 
have to propagate to a larger radius in order to overcome 
the effects of stretch. The minimum flame velocity and 
the flame temperature decrease with increasing dilution. 
Therefore, by increasing the dilution, one can obtain a 
mixture which will be quenched by endothermic reactions 
before it can overcome the effects of stretch. This ap- 
pears to be the case with the Hj : Oj : Nj/2: 1 : 13 mixture, 
which does not support spherical flame propagation. 

It has been shown experimentally that the addition of 
inert gases to confined fuel-air mixtures can cause the 
flammability limits to become narrower That is, by 
adding sufficient amounts of an inert gas, we can obtain 
a mixture that does not support flame propagation. Our 
calculations support this conclusion even though we do not 
include any heat or radical loss to the confining walls. 

The actual limits obtained might, however, depend on pa- 
rameters such as the method of initiation, the duration 
of energy deposition, and the geometry of the system. By 
systematically varying these parameters as well as the 
stoichiometry, the numerical simulations can be used to 
help answer questions regarding the existence of fundamen- 
tal flammability limits. 
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TIME-DEPENDENT SIMULATIONS OF LAMINAR FLAMES 
IN HYDROGEN-AIR MIXTURES 

1. Introduction 

In this paper, we present results of detailed numerical simulations of laminar flames in 
hydrogen-air mixtures. We use these simulations to discuss two major problems: (1) the 
relative importance of various diffusive transport processes on flame propagation in fuel- 
lean, stoichiometric, and fuel-rich mixtures, and (2) the time-dependent behavior of flames 
near the rich flammability limit. 

Various numerical simulation techniques have been applied to study a wide variety of 
problems involving the ignition, propagation and quenching of hydrogen flames. However, 
there have not been many attempts to use the detailed modeling approach to study quantita- 
tively the relative importance of various physical mechanisms involved in flame propagation. 
In one such effort [1], numerical simulations were used to evaluate the relative importance 
of thermal conduction and diffusion of species in specific methane and methanol flames. 
The calculations seemed to indicate that flame propagation is not governed primarily by 
diffusion of radicals but diffusion of fuel which subsequently affects the chemical reactions. 
More recently, we (2) showed that neither diffusion of species nor thermal conduction alone 
could give the correct burning velocity for a stoichiometric hydrogen-oxygen- nitrogen mix- 
ture. In the first part of this paper, we present a systematic study of the relative importance 
of various diffusive transport processes to flame propagation in fuel-lean, stoichiometric and 
fuel-rich hydrogen-air mixtures. 

In an earlier study on the effects of curvature and dilution on flame propagation [3]. 
we showed that a flame can be extinguished with less dilution in one geometry than in 
another, even in the absence of external heat losses and buoyancy effects. Specifically, 
we showed that a planar flame can propagate steadily in a dilute mixture which does not 
support a spherically expanding flame. Because this behavior is related to the effects of 
flame stretch, it raises the fundamental question of whether there is an extinguishment limit 
in the absence of stretch effects. This question can be addressed by carefully studying the 
behavior of planar flames in the absence of external heat losses and gravitational effects. 
Both lean and dilute hydrogen-oxygen-nitrogen mixtures are subject to flame instabilities. 
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Therefore, we restrict ourselves to the extinguishment behavior of fuel-rich hydrogen-air 
flames in the second part of this paper. 

2. The Numerical Model 

The results presented in this paper are derived from numerical simulations using 
FLAME1D, a one-dimensional, Lagrangian model [4,5]. The model solves the time- 

dependent conservation equations for mass, momentum and energy [6,7]. The model has 
been used for a variety of flame studies, including calculations of burning velocities [5], 
minimum ignition energies [8] and quench volumes [8]. The model has a modular form 
and permits a wide variety of geometric, initial, boundary, and time varying energy input 
conditions. 

The convective transport is solved by ADINC, a Lagrangian algorithm which solves 
implicitly for the pressure [9]. Because the method is Lagrangian and there is no numer- 
ical diffusion due to convective transport, the method gives an accurate representation of 
material interfaces and allows steep gradients in species and temperature to develop and 
be maintained. An adaptive Lagrangian regridding algorithm refines the grid ahead of the 
flame-front so that the flame zone remains well resolved throughout the calculation. Ther- 
mal conduction, ordinary diffusion, and thermal diffusion processes are included for each 
chemical species. The chemical interactions are described by a set of nonlinear, coupled or- 
dinary differential equations which are solved using a fully vectorized version of the selected 
asymptotic integration method CHEMEQ [10,11]. The algorithms representing the various 
chemical and physical processes are integrated separately and then asymptotically coupled 
by timestep-splitting techniques [7]. 

For this study we have used the hydrogen-oxygen reaction scheme [12] which involves 
the eight reactive species Hj, Oj, H, 0, OH, H^O, H02, HjOo, and the diluent, N4. The 
thermochemical properties of the various species involved are taken from the JANAF ta- 
bles [13]. The chemical reaction mechanism [4,8,12], consisting of about fifty steps has been 
extensively tested and shown to give good results. Burks and Oran [12] showed that the 
results computed with this mechanism compared well with experimentally observed induc- 
tion times, second explosion limits and the temporal behavior of reactive species. Oran et 
al. [14] have shown that when coupled to a convective transport algorithm, the mechanism 
gives good results in the simulation of the conditions behind a reflected shock. The lami- 
nar burning velocities calculated using the mechanism are in agreement with experimental 
data [5]. 

In the calculations presented below, the model was configured with an open boundary 
at one end to simulate an unconfined system. All the calculations were performed in a planar 
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(Cartesian) geometry. The initial temperature and pressure of the mixtures considered were 
298 K and 1 atm, respectively. 

The calculations were initialized with a Gaussian temperature profile. The central 
temperature and width were chosen so that the added initial energy was above the minimum 
ignition energy for the mixture. That is, it was ensured that the initial conditions provided 
enough energy for the flame to evolve to a steady profile. 

For either thin or planar flames, the instantaneous normal burning velocity can be 
calculated from the flame velocity if we know the velocity of the unburnt gases ahead of 
the flame. As shown in Fig. 1, for planar flames the velocity of the unburnt gases ahead of 
the flame is constant. Hence the burning velocity can be unambiguously determined as the 
difference betweeen the flame velocity and the flow velocity ahead of the flame, 

kiurn = V/lama — V/luid • 

We have discussed this definition and its application to curved flames elsewhere [3]. 

3. Relative Effects of Various Diffusive Transport Processes 

A series of calculations have been performed to study the effects of the individual 
diffusive transport processes on flame propagation in hydrogen-air mixtures. The three 
different diffusive transport processes we consider are 1) thermal conduction, 2) ordinary- 
diffusion, and 3) thermal diffusion. In order to isolate the relative importance of each, we 
have performed calculations with different combinations of these processes turned on. We 
repeated the same sequence of calculations for three mixtures: 

Fuel-lean: H 2 :0 2 :N 2 / 0.65:1:3.76 
Stoichiometric: H 2 :0 2 :N 2 / 2:1:3.76 
Fuel-rich: H 2 :0 2 :N 2 / 15:1:3.76. 

Table 1 summarizes the calculated burning velocities from this series of calculations. For 
each of the three mixtures considered we show four different cases: 

(a) All three diffusive transport processes are on; 

(b) Thermal diffusion is off, and thermal conduction and ordinary diffusion are on; 

(c) Thermal conduction is on, and thermal diffusion and ordinary diffusion are off: 

(d) Thermal conduction off, and ordinary and thermal diffusion are on. 

Several interesting trends in burning velocities appear from looking at this table. F'usr. 
we see that the relative importance of thermal diffusion does not change from lean to rich 
mixtures. However, the relative importance of thermal conduction and ordinary diffusion 
are different for lean, stoichiometric, and rich mixtures. 
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3.1 The Effect of Therma l Diffusion on Burninz_¥elocities 


By comparing cases (a) and (b) in Table 1, we see that the effect of thermal diffusion 
is to lower the burning velocities by 6-10% in all three mixtures. The effects of thermal 
diffusion in hydrogen-air mixtures have been studied earlier (4,15,16). 

Hydrogen atoms are the fastest moving, most easily diffused species in the system. 
Their diffusion behavior dominates the diffusion properties of the hydrogen-air system. 
Hydrogen atoms are formed in the high-temperature region of the flamefront, and are then 
diffused into the cold mixture by ordinary diffusion. Earlier studies have shown that thermal 
diffusion tends to diffuse hydrogen atoms from the lower temperature region to the higher 
temperature region [4,15]. Thus these two processes have the opposite effect with respect 
to moving hydrogen atoms. The increase in burning velocity when there is no thermal 
diffusion occurs because thermal diffusion inhibits the diffusion of hydrogen radicals into 
the cold, unburnt mixture ahead of the flamefront. 

3.2 Stoichio metric Mixture 

For the stoichiometric mixture, thermal conduction and ordinary diffusion are almost 
equally important in determining the burning velocity of the flame. Without thermal con- 
duction, the burning velocity is reduced to 62% of its standard value, i.e., the value of the 
burning velocity with all of the diffusive transport processes turned on. Without ordinary 
diffusion, it is reduced to 46% of its standard value. Figure 2 shows the behavior of the 
dame and burning velocities as a function of time. The initial condition causes a large initial 
velocity which eventually relaxes to a steady flame velocity. 

3.3 Fuel-Lean Mixture 

For the lean mixture, Table 1 shows that turning off either thermal conduction or ordi- 
nary diffusion lowers the burning velocity. In particular, we note that thermal conduction 
is much more important than ordinary diffusion in determining whether or not the flame 
propagates. When the thermal conduction is turned off, the flame does not reach a steady 
burning velocity at all. 

Figure 3 shows the calculated flame velocity and fluid velocity as a function of time 
for case (d), where the flame appears to die. The flame velocity decreases with time, and 
becomes practically identical to the fluid velocity. The result is a negligible burning velocity. 
Using a time-dependent model to study burning velocities allows us to see the flame die 
out in time. No changes in the initial conditions which we tried could keep this flame from 
eventually “dying out” in the wav shown in the figure. 
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For this lean flame, radicals can diffuse ahead of the flame front, but without thermal 
conduction, the temperature ahead of the flame never becomes high enough to sustain 
the reaction. The fuel deficit also means a deficit of hydrogen radicals, which hinders the 
reaction process. Thermal conduction alone, however, carries forward enough heat to create 
new radicals and sustain the process. Figure 4 shows the difference in the hydrogen radical 
profiles for cases (a), (c) and (d) at typical times in the calculation. 

3.4 Fuel- Rich Mixture 

For the fuel-rich mixture, ordinary diffusion is the most important process in deter- 
mining the burning velocity. In contrast to the fuel-lean flame, the fuel-rich flame does 
not propagate without ordinary diffusion. In the fuel-rich case, not enough heat can be 
transferred by conduction alone to sustain the reactions. However, because of the excess of 
hydrogen radicals, ordinary diffusion alone can spread enough radicals ahead of the front 
to initiate chemical reactions and sustain the flame. 

4. Flames in Fuel-rich Mixtures 

In the second set of calculations we restricted ourselves to flames in fuel-rich hydrogen- 
air mixtures. In these calculations all the diffusive transport processes discussed above were 
included. The amount of hydrogen in the mixture was systematically increased and with 
each concentration of hydrogen, the calculations were carried out untill a steadily propa- 
gating flame (if any) was observed. The burning velocities were calculated as discusssed 
earlier. The burning velocity and the temperature of the burned gases (flame temperature) 
for the various mixtures are summarized in Table 2. As expected, the burning velocity and 
the flame temperature decrease with increasing fuel concentration. However, a couple of 
interesting observations can be made from the data presented in the table. 

First of all, for a 80.8% hydrogen-air mixture (Hj^.'N? / 20:1:3.76 ), we obtain a 
steady burning velocity of about 8 cm/s. This mixture is beyond the experimentally ob- 
served rich- flammability limit under normal gravity conditions. This difference between the 
numerical and experimental observations may be due to the effects of gravity, heat losses 
to the confining walls in the experiments or multidimensional effects. As mentioned earlier, 
these effects have not been included in these numerical simulations. For example, we know 
from previous calculations [3] that curvature (or stretch) can reduce the burning velocity. 

The second interesting observation from the numerical simulations is that a 82.2% 
hydrogen-air mixture (H2:02:N2 / 22:1:3.76 ) does not support a steadily propagating flame. 
We tried different ignition sources but in all cases the burning velocity decreased to zero. 
The behavior of the transient flame was qualitatively similar to that shown in Fig. 3 for a 
dying flame. The simulations indicate that there is a flammability limit even in the absence 
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of stretch, buoyancy and external heat losses. Furthermore this limit is higher than that 
observed in the standard flammability limit tube under normal gravity conditions. 

A possible reason for this flammability limit can be seen in the trend of flame tempera- 
tures shown in Table 2. The flame temperature for a 20:1:3.76 mixture is just 970 K. With 
increase in fuel concentration the temperature would decrease further. At these tempera- 
tures, it is possible that the H atoms are consumed faster in chain terminating reactions 
than produced in chain branching reactions. This would result in a depletion of H atoms 
and a decrease in temperature. This could eventually lead to a situation where the energy 
released in exothermic reactions is not sufficient to balance the energy diffused by the trans- 
port processes and hence a steady flame is not possible. This possibility is currently being 
investigated. 

5. Summary and Conclusions 

We have examined two fundamental problems in premixed laminar flames using a 
detailed one-dimensional model. The model is time- dependent and can calculate the evo- 
lution to a steady-state, propagating flame from an initial energy source. Because it is 
time-dependent, the absence of a steady-state can also be revealed by it. 

In the first problem, we examined the relative importance of thermal conduction, 
thermal diffusion, and ordinary diffusion in determining the burning velocities of laminar 
hydrogen-air flames. Three cases were examined in detail: a fuel-lean mixture, a stoichio- 
metric mixture, and a fuel-rich mixture. For the fuel-lean and fuel-rich mixtures, care was 
taken to ensure that the mixture considered was far enough from the limit to have a finite, 
but relatively low, burning velocity. The results are summarized in Table 1, which shows 
the burning velocities calculated for the three mixtures. 

The effect of thermal diffusion is to decrease the burning velocity by 5-10%. in agree- 
ment with previous calculations [4,15,16]. This can be explained by noting that the effect 
of thermal diffusion is to diffuse hydrogen radicals in towards the hot region of the flame, 
instead of out towards the unburned region. It has the opposite effect of ordinary diffusion 
and thus inhibits the flame- propagation process. 

For the stoichiometric mixture, both thermal conduction and ordinary diffusion make 
important contributions to sustaining the flame. When thermal conduction was turned otf. 
the burning velocity was reduced to 62 % of its standard value. When ordinary diffusion 
was turned off, the burning velocity decreased to 46% of its standard value. These results 
were qualitatively as expected. 

The more surprising results were for the fuel-rich and the fuel-lean mixtures. For tin' 
fuel-lean mixture, the flame would not propagate when thermal conduction was turned oiF. 
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That is, ordinary diffusion alone could not sustain the flame. However, the burning velocity 
was only 69% of its standard value with just thermal conduction. 

In contrast to the fuel-lean behavior, the fuel-rich mixture could not sustain a flame 
when ordinary diffusion was turned off. That is, thermal conduction alone could not sustain 
the flame. With just ordinary and thermal diffusion, the burning velocity was 86% of its 
standard value. 

Our general conclusion from this study is that both ordinary diffusion and thermal 
conduction are necessary to quantitatively describe flame propagation in a hydrogen-air 
mixture. Their relative importance, however, varies as we go from fuel-lean to fuel-rich 
hydrogen-air mixtures. It is not clear that these conclusions hold for hydrocarbon combus- 
tion, since there are other, heavier radicals that might change the balance between thermal 
conduction and ordinary diffusion. An implication of this work to flame modelling efforts 
is that it may be possible to neglect thermal diffusion, but both molecular diffusion and 
thermal conduction must be included in the model. 

In the second part of this paper we considered the behavior of flames in fuel rich 
hydrogen-air mixtures near the experimentally observed flammability limit. The effects of 
gravity, stretch and external heat losses were eliminated in the numerical simulations. The 
results suggest wider flammability limits than those observed experimentally under normal 
gravity conditions. The simulations indicate that there may be a limit due to chemical 
kinetic considerations alone. Further details of a possible chemical-kinetic limit are currently 
under investigation. The simulations also suggest that multi-dimensional calculations under 
normal gravity conditions and experiments under reduced gravity conditions are needed to 
better understand flames near the standard flammability limits. 
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Table 1. - Burning Velocities of Hj-Air Mixtures 


Case 

Processes On 

Lean 

Stoichiometric 

Rich 

(a) 

All 

36 cm/s 

185 cm/s 

44 cm/s 

(b) 

Thermal Conduction 
Ordinary Diffusion 

40 

195 

48 

(c) 

Thermal Conduction 

25 

85 

0 

w 

Ordinary Diffusion 
Thermal Diffusion 

0 

115 

38 


Table 2. - Flames in Fuel* Rich Mixtures 


Mixture 
Hi : Oi : /Vo 

Burning Velocity 
(cm/s) 

Flame Temperature 

(K) 

16:1:3.76 

38 

1130 

17:1:3.76 

30 

1108 

18:1:3.76 

23 

1046 

19:1:3.76 

16 

1002 

20:1:3.76 

8 

970 

22:1:3.76 

— 

<970 
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POSITION (mm) 

Velocity and temperature profiles in a planar flames propagating in an H 2 :0 2 :N 2 /2:I4 mixture. 




VELOCITY (cm/s) 



TIME (ms) 

Fig. 2 — Calculated fluid, flame, and burning velocities as a function of time 
for the mixture H2:0 3 N 2 /2:1:3.76 at 298 K and 1 atm. 
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Fig. 3 — Calculated flame and fluid velocities as a function of time 
for the H 2 :O 2 :N 2 /0.65: 1:3.76 mixture. 
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Abstract 


Time-dependent two-dimensional simulations of perturbed premixed laminar flames have 
been used to study the development of cellular structures in rich and lean hydrogen flames. The 
model includes a detailed reaction mechanism for hydrogen-oxygen combustion consisting of 25 
elementary reactions among eight reactive species and a nitrogen diluent, molecular diffusion 
between all of the species, thermal conduction, and convection. Calculations of the evolution 
of a perturbed lean hydrogen flame showed that the flame was unstable at the front, and 
the structure that evolved resembled the cellular structures observed in experiments. The same 
perturbation applied to a rich hydrogen flame showed that the perturbation died out and cellular 
structures did not appear. Tests of the importance of molecular diffusion were carried out by 
varying the binary diffusion coefficients. First, the binary diffusion coefficient of molecular 
hydrogen was set equal to that of molecular oxygen, in which case the perturbation died out. 
Then, the binary diffusion coefficient of molecular oxygen was set equal to that of molecular 
hydrogen, in which case the instability persisted. The results of the simulations support the 
diffusional-thermal theory of flame instability. 


I. Introduction 

Previous one-dimensional, time-dependent numerical simulations of laminar premixed 
flames have given physical insights and quantitative information on the effects of parame- 
ters such as stretch, curvature, and dilution on flame initiation and propagation, as well as 
qualitative information on flammability limits 1,3 . The structure of flames, especially near the 
flammability limits, is usually multidimensional. To investigate such multidimensional effects, 
we have developed a detailed time-dependent two-dimensional model that can be used to sim- 
ulate either premixed or diffusion flames. In this paper, we present results of calculations that 
used this model to study premixed laminar flames in rich and lean hydrogen-oxygen mixtures 
diluted with nitrogen. The results of these simulations are discussed in terms of current theo- 
ries of flame stability and the formation of cellular structures. The validity of two prominent 
theories of cellular instability is investigated. 

Early experimental observations showed that cellular flames occur in lean hydrogen-air 
as well as in rich hydrocarbon-air mixtures such as propane and ethane 3,4 . Cellular flames 
have generally been observed in lean fuel-air mixtures when the fuel is lighter than air and 
in rich fuel-air mixtures when the fuel is heavier than air 4,s . These observations suggest that 
the relative diffusivities of fuel and air play a crucial role in determining the sensitivity of a 
mixture to cellular instability and lead to the formulation of the preferential diffusion theory 5 . 
According to this theory, preferential diffusion of the lighter reactant will create regions of 
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varying stoichiometry in front of a non-planar flame. This in turn causes portions of the flame 
to move at different speeds, making the flame front either less or more planar. This theory, 
however, does not explain all of the observations of cellular flames. For example, ethane and 
ethylene have molecular weights close to but less than that of oxygen, and yet they show cellular 
structure in rich mixtures. 

Another theory of cellular flames is the diffusional- 1 hermal theory 6,7 . In this theory, an 
abundance of the excess component was assumed, so that the reaction is controlled only by the 
deficient component. Then for the simple case of a one-step reaction, in what is effectively a 
single reactant system, this theory predicts that cellular instability occurs when the thermal 
diffusivity of the mixture is sufficiently small compared to the mass diffusivity of the reactant. 
For lean hydrogen-air mixtures, hydrogen is the limiting reactant and its mass diffusivity sig- 
nificantly exceeds the thermal diffusivity of the mixture. In rich hydrogen-air mixtures, oxygen 
is the limiting reactant and its mass diffusivity is nearly the same as the mixture thermal diffu- 
sivity. Hence the theory agrees with the early experimental predictions that lean hydrogen-air 
mixtures are unstable and rich mixtures are stable. 

A limitation of the diffusional-thermal theory is that it assumes that the density of the fluid 
is a constant, and hence the effects of thermal expansion are neglected. Furthermore, according 
to this theory, the cell size increases indefinitely as the Lewis number goes to its critical value. 
More recent experimental observations 6 * 9 have shown cellular instabilities in some rich and near 
stoichiometric hydrogen-air mixtures, in contrast to the predictions of both the preferential 
diffusion as well as the diffusional-thermal theories. The diffusional-thermal theory is strictly 
valid only for strongly non-stoichiometric mixtures, although it has been extended 10,11 to near 
stoichiometric mixtures by considering both the deficient and abundant components. According 
to this two-reactant theory, the stability behavior depends on the stoichiometry of the mixture, 
the difference between the Lewis numbers of the two components, and the reaction orders of 
each reactant. This observation suggests that in multicomponent systems, the Lewis numbers 
of many of the components and detailed chemical kinetics might play a role. The modified 
diffusional-thermal theory also neglects the thermal expansion of the combustion products. 

The simulations of multidimensional flame structure presented here include a multi-reaction 
mechanism for hydrogen combustion, molecular diffusion between the reactants, intermediates, 
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and products, thermal conduction, and convection. Such a detailed model should allow us to 
investigate the tendency of mixtures to show cellular instability, to evaluate the importance of 
various contributing physical processes, and to evaluate the predictions of the various theories. 
For example, previous explanations of cellular instability have considered the fuel and oxidizer, 
but the instability may also depend on the diffusivities of intermediate radicals. Lh addition, 
the effects of the fluid dynamics and therefore the variable gas density can be evaluated. Below 
we describe the numerical model, show the results of using this model to simulate a series 
of hy d r ogen-oxy gen flames with different stoichiometries, and then relate these results to the 
predictions of two theories of flame instability. 

II. Multidimensional Flame Model 

A detailed model of a flame must contain accurate representations of the convective, dif- 
fusive, and chemical processes. The individual importance of these processes varies from rich 
to lean flames, and is especially notable near the flammability limits 3 where the exact behavior 
of these flames depends on a delicate balance among the processes. In this section we briefly 
describe the algorithms and input data used to model and couple the effects. 

The fluid convection algorithm must be able to maintain the sharp gradients present in 
flames. Numerically this means that any important numerical diffusion in the calculation must 
be considerably less than any physical diffusion effect. Many explicit algorithms now exist that 
treat sharp discontinuities in flow variables accurately, but these methods are extremely ineffi- 
cient at the very low velocities associated with laminar flames. The Barely Implicit Correction 
Flux- Corrected Transport (BIC-FCT) algorithm 13 was developed specifically to solve low-speed 
flow problems with high accuracy. BIC-FCT solves the coupled continuity equations for a vari- 
able density flow. BIC-FCT combines an explicit high-order, nonlinear FCT method 13,14 with 
an implicit correction process. This combination maintains high-order accuracy and yet re- 
moves the timestep limit imposed by the speed of sound. By using FCT for the explicit step, 
BIC-FCT is accurate enough to compute with sharp gradients without overshoots and under- 
shoots. Thus spurious numerical oscillations that would lead to unphysical chemical reactions 
do not occur. The computational cost per timestep needed for BIC-FCT is about the same as 
for the standard explicit FCT method 13 , which is an amazing gain considering that BIC-FCT 
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is implicit and the standard FCT is explicit. 

Thermal conductivity of the individual species is modeled by a polynomial fit in tempera- 
ture to existing experimental data. Individual conductivities are then averaged using a mixture 
rule 1516 to get the thermal conductivity coefficient of the gas mixture. A similar process is 
used to obtain the mixture viscosity from individual viscosities. Heat and momentum diffusion 
are then calculated explicitly using these coefficients. In the problem considered in this paper, 
the timestep restriction of an explicit method for the diffusion terms does not cause any loss in 
efficiency. 

Mass diffusion also plays a major role in determining the properties of laminar flames. 
Binary mass diffusion coefficients are represented by an exponential fit to experimental data, 
and the individual species diffusion coefficients are obtained by applying mixture rules 19 . The 
individual species diffusion velocities are solved for explicitly by applying Fick’s law followed 
by a correction procedure to ensure zero net flux 19 . This procedure is equivalent to using the 
iterative algorithm DFLUX 17 which solves the multi-component diffusion equations exactly. 
This method is substantially faster than one that uses matrix inversions and is well suited for a 
vector computer. This algorithm is also explicit, but because the effective Lewis number of the 
mixture is close to one, the timestep suitable for heat conduction is adequate for mass diffusion 
as well. 

Chemistry of the hydrogen-oxygen flame is modeled by a set of 24 reversible reactions 
describing the interaction of eight reacting species, Hj, Oj, H, 0, OH, HOa, HjOj, HjO, and 
Nj as a non- reacting diluent 18 . This reaction set is solved at each timestep with a vectorized 
version of CHEMEQ, an integrator for stiff ordinary differential equations 19 . Because of the 
complexity of the reaction scheme and the large number of cells in a two-dimensional calculation, 
the solution of the chemical rate equations takes a large fraction of the total computational 
time. A special version of CHEMEQ called TBA was developed to exploit the special hardware 
features of the CRAY X-MP vector computer. The full potential of new computer architectures 
can only be tapped at the expense of portability of programs to other machines. 

All of the chemical and physical processes are solved sequentially and then are coupled 
asymptotically by timestep splitting 30 . This modular approach greatly simplifies the model 
and makes it easier to test and change the model. Individual modules were tested against 
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known analytic and other previously verified numerical solutions. One- dimensional predictions 
of the complete model were compared to those from the Lagrangian model FLAME1D which 
has been benchmarked extensively against theory and experiment 15 . The exact values of the 
input data for a hydrogen-oxygen flame, for example, the chemical reaction rates, thermody- 
namic parameters, and molecular diffusion coefficients, have been described in detail in previous 
papers 15 , and therefore for brevity are only referenced here. 

III. Tests of Mechanisms of Laminar Flame Instability 

Initial conditions for the two-dimensional calculations were taken from one-dimensional 
calculations run to give the conditions for steady, propagating flames. Figure 1 describes the 
configuration under study and gives the boundary conditions of the computational domain. 
Fresh unburnt gas flows in from the left, and the products of chemical reaction at the flame 
front flow out at the right. If the inlet velocity is set to the burning velocity of the flame, 
the flame zone is fixed in space and a steady solution is obtained. Thus, the transient effects 
arising from the ignition process can be eliminated and the one-dimensional solution provides 
the initial condition for the two-dimensional calculation. The two-dimensional computational 
domain was 2.0 cm x 4.5 cm, which was resolved by a 56 X 96 variably spaced grid. Fine zones, 
0.36 mm x 0.15 mm, were clustered around the flame front. 

Flames in a Hydrogen-Lean Mixture 

The first calculation is of a flame in a fuel-lean mixture of hydrogen and oxygen diluted with 
nitrogen, Ha:Oj:Nj / 1.5:1:10, a flame that was clearly multidimensional in the experiments by 
Mitani and Williams 9 . The initial condition described by Fig. 1 is perturbed by displacing the 
center portion of the flame in the direction of the flow. The frames in Fig. 2 show isotherms just 
after the perturbation and then their evolution at subsequent times. The isotherms indicate 
that the temperature increases in the center portion of the flame which is convex to the flow and 
decreases in the the two adjacent concave regions. This indicates that the reaction progresses 
more vigorously in the convex region. This conclusion is corroborated by the atomic hydrogen 
number density contours in Fig. 3. The concentration of atomic hydrogeii increases in the 
convex region and decreases in the concave regions. Also, the burning velocity in the convex 
region appears to be slightly higher than in a planar region, and in the concave regions the 
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burning velocity has noticeably decreased. We conclude that this lean one-dimensional dame is 
unstable, and a pattern resembling a cellular structure has appeared by the time the calculation 
was terminated. 

Flames in a Hydrogen-Rich Mixture 

The second case considered the evolution of a disturbance in a fuel-rich mixture of hydrogen 
and oxygen, HjcOj.Na / 3.0:1:16. This mixture is stable in the experiments of Mitani and 
Williams 9 , and also in our calculations. The initial disturbance is quickly damped, restoring 
the flame profile to its one-dimensional shape. The frames in Fig. 4 show the evolution of 
temperature contours for the rich flame. The perturbed center section rapidly straightens out, 
indicating that the flame speed here is lower than in the planar sections. Figure 5 gives the 
atomic hydrogen contours for the rich flame. The number density of atomic hydrogen is lower 
in the center, indicating a weaker reaction and hence a lower burning velocity. This is the 
opposite of the behavior we observed in the lean flame calculations. 

These two results are in agreement with experiments and with both the theory of pref- 
erential diffusion and the diffusional-thermal theory. In the next two studies, the diffusion 
coefficients of the reactants (Hj and Oj) were changed in a systematic manner in order to 
determine the correct role of mass diffusion. All other variables were held the same as in the 
fuel-lean case (HjiOjiNj / 1.5:1:10). 

The Role of Mass Diffusion 

By setting the diffusion coefficient of hydrogen equal to that of oxygen, the light reactant 
no longer moves faster than the heavier reactant. This change in the model resulted in a stable 
one-dimensional flame. Figure 6 shows that the isotherms straighten out and slowly return to 
the unperturbed planar condition. A lower concentration of atomic hydrogen was found in the 
center section of the flame front, which agrees with the observed stability of the flame. 

The fact that the modified lean flame described above is stable supports the theory of 
preferential diffusion. Because the diffusion coefficient of hydrogen equals that of oxygen, there 
can be no preferential diffusion of hydrogen to cause the instability. The remaining stabilizing 
factors, such as heat conduction, restore the flame to its unperturbed position. However, these 
results are also in agreement with the diffusional-thermal theory. Since the Lewis number of 
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the reactants in the modified lean flame is dose to unity, the diffusionai-thermal mechanism for 
instability is greatly reduced. The stabilizing factors overcome any tendency to instability. 

A further test was needed to distinguish between these two mechanisms. In this test, 
the diffusion coeffitients of the oxygen molecule were set to that of the hydrogen molecule. 
Now there are two fast species. All other conditions were unchanged from the lean flame test. 
This case again eliminates preferential diffusion and thus if preferential diffusion is indeed the 
mechanism, this flame should be stable. However, as Fig. 7 indicates, the flame is unstable 
and dosely resembles the standard lean-flame result shown in Fig. 4. This result contradicts 
the predictions of the preferential diffusion theory. The diffusionai-thermal theory, on the other 
hand, predicts that because the mass diffusion of the defident spedes, molecular hydrogen, 
remains high, this flame will be unstable. Thus the series of numerical simulations presented 
here support the diffusionai-thermal theory of cellular instability. 

IV. Condusion 

A new time-dependent two-dimensional model has been used to study the stability of rich 
and lean premixed laminar hydrogen-oxygen flames diluted with nitrogen. The model included 
a detailed chemical reaction scheme consisting of 24 reversible reactions among eight species, 
thermal conductivity, and the individual spedes diffusion velodties. Four tests were presented 
which showed the evolution of perturbed two-dimensional flames: 

1. A lean flame, H]:02:N2 / 1.5:1:10, 

2. A rich flame, H2:02:N2 / 3.0:1:16, 

3. The lean flame with the H 2 diffusion coefficients set equal to the O 2 diffusion coeffi- 
cients, and 

4. The lean flame with the O 2 diffusion coeffidents set equal to the H 2 diffusion coeffi- 
cients. 

The numerical simulations show that the one-dimensional lean flame is unstable, and by the 
time the calculation is terminated, evolves into a pattern that resembles a cellular flame. The 
rich flame, however, is stable, and the perturbation dies in time. Both of these results agree with 
those in the experiments of Mitani and Williams 9 for mixtures with the same stoichiometry. 

We then ask if we can use the simulation to determine the mechanism controlling the 
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instability. This is addressed by the third and fourth calculations, which focus on the mass 
diffusion of the reactants, Hj and 0 2 . In the third calculation, the mass diffusivity of hydrogen 
was set equal to oxygen. This eliminates the preferential diffusion effect, and means that the 
light fuel will not move faster than the other reactant. If the mechanism of instability were 
preferential diffusion, the lean flame would be stable. This is in fact the result: the lean flame 
in the calculation with this change of relative diffusivities was stable. 

This result supports both the preferential diffusion theory and the diffusional-thermal 
theory. A further calculation was performed in which the mass diffusivity of oxygen was set 
equal to that of hydrogen. In this case, the flame was unstable, which disagrees with the 
prediction of preferential diffusion. However, the diffusional-thermal theory does predict that 
this flame is unstable. Thus for these hydrogen flames studied, the results agree with the 
predictions of the diffusional-thermal theory proposed by Barenblatt* and Sivashinsky 7 . 

This series of calculations has also raised a number of computational issues. First, there 
were limitations on the calculations which it would be interesting to eliminate. We would in fact 
like the computational domain to be larger with the same or even better resolution. A larger 
domain would allow the central parts of the calculation to be less influenced by the boundary 
conditions, and it would be interesting to see if more cells form and how they change in time. 
It would be useful to redo these calculations with other kinds of perturbation at the front, for 
example, vary the wavelength of the perturbation. This would be used to address questions such 
as what is the necessary width of the perturbation relative to the flame thickness to cause the 
lean flame to go unstable, and what is the effect of long versus short wavelength perturbations 
on the growth and type of instability that results. 

The computations were very expensive. The cost of each one was between four or five hours 
of CRAY X-MP 2/4 time, which translates to 0.5 ms per pointstep. Because the computations 
were so expensive, calculations were terminated even though it would have been interesting to 
follow them longer. For example, we would like to address the question of whether the cellular 
structure in Fig. 2 will split into more cells. None of the computations were bound by the 
memory limitation of the computer, but rather were CPU bound, with most of the time spent 
in integrating the extensive chemical reaction mechanism. One approach we could take, now 
that we have been able to reproduce the types of results that are observed experimentally. 
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is to consider a model system with a greatly simplified chemical reaction mechanism. This 
appears justified for the study of the mechanism of cellular structure because chemistry does 
not appear to play a major role in the instability. Then the numerical model would be relatively 
inexpensive to run, and could be used to test selected theoretical concepts much less expensively 
and perhaps less ambiguously. 

There are, however, certain advantages to having the full simulation model. For example, 
other phenomena we wish to investigate, that occur near the flammability limits, might be very 
sensitive to the details of the chemical reactions. The mass diffusivities of light intermediate 
species such as atomic hydrogen might also play a role in determining stability. At this juncture, 
now that the model exists, it is possible to move in either direction. 
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Figure 1 . Initial and boundary conditions for the two-dimensional flame calculations. 
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Figure G. Involution of isotherms, dilfusioa coefficients of II 2 equal to 0; 
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